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dense  slope  representations  as  well  as  dense  depth  representations.  This 
marked  change  from  previous  surface  reconstruction  algorithms  allows 
a  natural  integration  of  orientation  constraints  into  the  surface  description, 
a  feature  not  easily  incorporated  into  earlier  algorithms.  In  addition, 
the  coupled  depth/slope  model  generalizes  to  allow  for  varying  amounts 
of  smoothness  at  different  locations  on  the  surface. 

This  computational  model  helps  conceptualize  the  problem  and  leads  to 
two  possible  implementations  -  analog  and  digital.  The  model  can  be 
implemented  as  an  electrical  or  biological  analog  network  since  the  only 
computations  required  at  each  locally  connected  node  are  averages,  additions 
and  subtractions.  A  parallel  digital  algorithm  can  be  derived  by  using 
finite  difference  approximations  The  resulting  system  of  coupled  equations 
can  be  solved  iteratively  on  a  m4sh-of -processors  computer,  such  as  the 
Connection  Machine.  Furthermore,  concurrent  multi-grid  methods  are 
designed  to  speed  the  convergence  of  this  digital  algorithm. 
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Abstract 


Reconstructing  a  surface  from  sparse  sensory  data  is  a  well  known  prob¬ 
lem  in  computer  vision.  Early  vision  modules  typically  supply  sparse  depth, 
orientation  and  discontinuity  information.  The  surface  reconstruction  module 
incorporates  these  sparse  and  possibly  conflicting  measurements  of  a  surface  into 
a  consistent,  dense  depth  map.  The  coupled  depth/slope  model  developed  here 
provides  a  novel  computational  solution  to  the  surface  reconstruction  problem. 
This  method  explicitly  computes  dense  slope  representations  as  well  as  dense 
depth  representations.  This  marked  change  from  previous  surface  reconstruc¬ 
tion  algorithms  allows  a  natural  integration  of  orientation  constraints  into  the 
surface  description,  a  feature  not  easily  incorporated  into  earlier  algorithms.  In 
addition,  the  coupled  depth/siope  model  generalizes  to  allow  for  varying  amounts 
of  smoothness  at  different  locations  on  the  surface. 

This  computational  model  helps  conceptualize  the  problem  and  leads  to  two 
possible  implementations  -  analog  and  digital.  The  model  can  be  implemented 
as  an  electrical  or  biological  analog  network  since  the  only  computations  re¬ 
quired  at  each  locally  connected  node  are  averages,  additions  and  subtractions. 
A  parallel  digital  algorithm  can  be  derived  by  using  finite  difference  approxima¬ 
tions.  The  resulting  system  of  coupled  equations  can  be  solved  iteratively  on 
a  mesh-of-processors  computer,  such  as  the  Connection  Machine.  Furthermore, 
concurrent  multi-grid  methods  are  designed  to  speed  the  convergence  of  this 
digital  algorithm. 
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CHAPTER  1 


INTRODUCTION 


This  report  studies  surface  reconstruction,  a  well  known  problem  in  com¬ 
puter  vision.  Computer  vision  is  an  active  topic  of  current  research  for  two 
reasons.  First,  we  want  to  build  artificial  vision  systems  for  use  in  robotics,  fac¬ 
tory  automation,  parts  inspection  and  other  industrial  applications.  Second,  we 
study  computer  vision  in  order  to  understand  the  information  processing  tasks 
necessary  for  biological  vision  systems. 

This  computational  approach  to  vision  was  first  proposed  by  Marr  [  1982] . 
Marr  partitioned  visual  processing  into  three  levels:  early  vision,  in  which  various 
modules  detect  shape,  depth  or  discontinuities  from  raw  images;  intermediate 
vision,  which  uses  the  outputs  from  early  vision  modules  and  reconstructs  sur¬ 
faces  in  a  viewer-centered  coordinate  frame,  and  late  vision,  in  which  processes 
such  as  recognition  are  performed  on  objects  in  an  object-centered  coordinate 
frame. 

In  the  years  since  Marr’s  original  work,  a  tremendous  amount  of  research  has 
been  done  in  early  vision.  Underlying  much  of  this  work  is  the  concept  that  early 
vision  algorithms  exhibit  extraordinarily  fine-grained  locality  and  thus  fit  neatly 
onto  the  coming  wave  of  massively  parallel  computers.  It  has  not  been  clear, 
however,  how  the  speedups  achievable  in  early  vision  algorithms,  through  the 
use  of  parallel  computation,  can  be  extended  to  later  levels  of  visual  processing. 

This  report  concentrates  on  increasing  the  processing  speed  for  intermediate 
vision  by  developing  a  novel  way  of  conceptualizing  surface  reconstruction.  This 
new  approach,  called  the  coupled  depth/slope  model,  lends  itself  to  implementa- 
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an  analog  network  realization  for  an  even  faster  solution.  Analog  network  ideas 
are  not  new  to  computer  vision.  Horn  [1974]  designed  an  analog  network  for 
the  solution  of  the  lightness  from  brightness  problem.  Poggio  and  Koch  [1985] 
used  regularization  analysis  to  devise  analog  hardware  to  solve  most  of  the  prob¬ 
lems  in  early  vision.  Here  we  develop  an  analog  network  model  to  solve  the  full 
surface  reconstruction  problem  in  intermediate  vision. 

Perhaps  an  even  more  interesting  result  of  this  work  is  the  emerging  resem¬ 
blance  of  the  artificial  vision  systems  we  are  constructing  to  the  biological  vision 
systems  that  we  are  studying.  This  convergence  of  the  two  research  directions, 
although  not  required,  is  not  coincidental  either. 

Necessary  criteria  for  biologically  feasible  information  processing  systems 
have  been  compiled  by  [Ullman  1979,  p.  88;  Grimson  1981,  p.  163;  Knight 
1983,  p.  7]  and  can  be  summarized  as  follows: 

1.  Massive  parallelism  -  a  large  number  of  independent  nodes  performing  com¬ 
putations. 

2.  Uniform,  simple  nodes. 


3.  Local  interconnect  scheme  between  nodes. 

4.  Fault  tolerance  -  performance  does  not  depend  critically  on  any  one  node. 

These  same  constraints  for  biological  feasibility  carry  over  to  the  domains 
of  VLSI  design  and  the  construction  and  programming  of  massively  parallel  ar¬ 
chitectures.  In  VLSI  design,  parallelism  on  chip  provides  more  efficient  use  of 
real  estate;  uniform,  simple  cells  shorten  layout  time;  local  interconnect  reduces 
the  need  for  long  wires,  and  fault  tolerance  increases  yield,  as  a  chip  can  still 
function  even  though  some  parts  are  broken.  In  parallel  architecture  design  and 
use.  highly  parallel  algorithms  keep  the  maximum  number  of  processors  busy; 
uniformity  in  both  hardware  and  software  shortens  the  time  for  design,  testing 


and  debugging;  local  communication  in  algorithms  allows  for  simpler  communi¬ 
cation  hardware,  and  fault  tolerance  enables  the  machine  to  keep  running  even 
though  some  processors  or  connections  may  become  disabled.  These  same  con¬ 
straints  that  hamper  VLSI  and  computer  architecture  designers  seem  to  exist 
in  the  biological  world  as  well.  Since  these  three  implementation  mediums  are 
bound  by  similar  constraints,  it  should  not  be  surprising  that  similar  algorithms 
result. 

Parallelism,  locality,  simplicity  and  uniformity  are  features  of  the  algorithms 
developed  in  this  report.  This  work  concentrates  on  the  artificial  machine  vision 
approach  to  computer  vision.  We  try  to  reconstruct  the  original  surface  as 
fast  and  as  accurately  as  possible  in  both  digital  and  analog  hardware.  Some 
psychophysical  tests  have  shown  that  there  exist  some  similarities  between  the 
coupled  depth/slope  approach  and  human  vision.  These  similarities  are  welcome 
and  noted,  but  are  not  required. 

Through  surface  reconstruction  we  seek  to  find  a  plausible  surface  given 
sparse,  noisy  and  possibly  conflicting  measurements  of  an  original  surface  [Grim- 
son  1981;  Terzopoulos  1984].  Surface  interpolation  and  surface  approximation 
are  two  special  cases  of  surface  reconstruction.  A  surface  interpolation  algorithm 
fits  a  smooth  surface  over  sparse  but  exact  depth  measurements.  Surface  approx¬ 
imation  extends  this  idea  to  fit  surfaces  over  noisy  depth  measurements.  The 
final  computed  surface  is  not  required  to  pass  through  the  data  measurements. 
Surface  reconstruction  is  a  more  general  problem,  combining  information  from 
many  modalities.  This  report  considers  the  integration  of  depth  and  orientation 
constraints  with  surface  discontinuity  information.  A  solution  of  this  problem  is 
extremely  computationally  intensive.  Chapter  2  describes  surface  reconstruction 
more  fully. 

Chapter  3  develops  the  new  coupled  depth /slope  model  for  surface  recon- 
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struction.  This  network  model  is  shown  to  be  well  suited  for  the  surface  recon¬ 
struction  problem.  The  coupled  depth/slope  model  is  an  improvement  on  the 
thin  plate  model  used  by  Terzopoulos  [1984]  and  other  biharmonic  models  for 
the  following  reasons: 

1.  Orientation  constraints  are  integrated  into  the  surface  solution  very  nat¬ 
urally.  Other  models,  such  as  the  thin  plate  model,  contain  no  inherent 
notion  of  slope. 

2.  The  model  generalizes  to  arbitrary  levels  of  smoothness.  In  fact  it  directly 
implements  the  controlled  continuity  models  of  Terzopoulos  [1986]. 

3.  Both  digital  and  analog  network  implementations  of  the  coupled  depth /slope 
model  are  possible  and  both  are  discussed  in  this  report.  The  parallel  ver¬ 
sion  of  the  standard  biharmonic  mask  used  by  Terzopoulos  and  Grimson 
does  not  converge;  an  extra  local  weighting  step  must  be  added  to  regain 
convergence.  The  parallel  version  of  the  coupled  depth/slope  algorithm 
described  here,  converges  by  simply  alternating  between  depth  and  slope 
calculations. 

4.  Both  the  analog  and  digital  implementations  of  this  model  require  only  near¬ 
est  neighbor  communication  throughout  the  entire  algorithm.  This  com¬ 
munication  feature,  plus  the  fact  that  the  nodes  are  performing  extremely 
simple  arithmetic  and  averaging  operations,  simplifies  both  the  analog  and 
the  parallel  digital  implementations. 

Chapter  4  derives  the  finite  difference  equations  and  the  resulting  iterative 
scheme  necessary  for  a  digital  algorithm.  This  algorithm  is  well  suited  for  a 
fine-grain  massively  parallel  computer,  such  as  the  Connection  Machine  [Hillis 
1985).  The  single-grid  digital  implementation  of  this  model  suffers  from  slow 
convergence  rate  problems  common  to  all  local  iterative  methods.  We  propose 
a  simple  but  fast  extension  to  analog  network  theory  to  speed  up  convergence 


through  the  use  of  coarser  levels  of  resolution.  This  idea  is  inspired  by  the 
concurrent  multigrid  work  of  Kuo  [1985]  and  Terzopoulos  [1984],  but  is  a  much 
simpler  realization.  This  multigrid  network  is  useful  in  speeding  up  both  the 
analog  and  digital  realizations  of  the  coupled  depth/slope  model. 

Ideas  for  an  analog  implementation  of  the  coupled  slope/depth  network  are 
detailed  in  Chapter  5.  Analog  networks  are  interesting  for  a  number  of  rea¬ 
sons.  Most  important  for  this  work  is  that  these  networks  are  extremely  fast. 
For  example,  they  can  compute  solutions  to  the  surface  reconstruction  prob¬ 
lem  several  orders  of  magnitude  faster  than  the  corresponding  parallel  digital 
algorithms.  Also  analog  networks  are  interesting  because  the  algorithms  they 
implement  may  be  similar  to  the  ones  used  by  biological  systems  for  vision. 
This  raises  the  possibility  for  someday  “growing”  and  “programming”  neurons 
to  create  special  purpose  vision  systems.  Chapter  6  discusses  promising  research 
directions  based  upon  this  work. 
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CHAPTER  2 


BACKGROUND 


Surface  reconstruction  is  an  integral  part  of  the  approach  to  vision  described 


by  Marr  [1982].  Marr  proposed  several  intermediate  stages  in  his  representa¬ 
tional  framework  for  deriving  shape  information  from  images.  His  representa¬ 


tional  stages  are  as  follows: 


Greyscale  images  -  arrays  of  intensity  information  presented  as  the  input  to 


the  visual  system. 


Edges  -  the  primal  sketch  makes  explicit  the  sharp  changes  of  intensity  in 


images. 


2|-D  sketch  -  an  explicit  representation  of  the  depth  and  orientation  at 


each  point  on  the  visible  surfaces,  plus  discontinuities  in  its  smoothness. 
The  2^-D  sketch  is  taken  from  the  same  viewer  centered  coordinate  system 
where  the  original  images  were  sampled. 


The  3-D  model  representation  -  an  object  centered  coordinate  system  that 


represents  explicitly  the  volumes  of  space  that  objects  occupy. 


Early  vision  modules  include  the  stereo  module,  the  “shape  from”  modules, 
and  the  edge  detection  module,  which  provide  sparse  and  possibly  conflicting 


sets  of  data  for  input  to  the  surface  reconstruction  module.  For  example,  stereo 
[Marr  and  Poggio  1979]  supplies  an  extremely  sparse  set  of  depth  values,  the 
shape  from  shading  module  (Horn  1975]  provides  local  surface  orientations,  and 
edge  detection  [Hildreth  1980;  Canny  1983]  provides  possible  locations  for  surface 
discontinuities.  The  relationship  of  the  surface  reconstruction  module  to  these 


earlier  stages  is  diagrammed  in  Figure  2-1. 


wasi 


In  tliis  report,  these  early  vision  stages  arc  considered  as  black  boxes.  We  do 
not  look  inside  them  to  sec  how  they  work;  we  only  process  their  outputs  with  the 
surface  reconstruction  module.  We  do,  however,  open  up  the  surface  reconstruc¬ 
tion  black  box  and  explain  its  operation.  In  the  future,  we  will  want  to  blur  the 
clean  distinction  between  early  and  intermediate  vision  to  improve  the  overall 
system.  For  example,  the  stereo  module  and  the  surface  reconstruction  module 
would  both  produce  more  accurate  results  if  there  were  feedback  between  the 
two  systems.  There  is  also  evidence  that  the  early  vision  modules  interact  with 
one  another  [Poggio  1984|.  The  discontinuity  detector  would  operate  much  bet¬ 
ter  if  it  had  a  closer  interface  to  the  stereo  module,  to  the  surface  reconstruction, 
and  to  other  visual  cues  such  as  structure  from  motion.  Nonetheless,  placing 
a  clear  dividing  line  between  early  vision  modules  and  surface  reconstruction 
sufGccs  for  this  report. 

Grimson  [1981]  implemented  the  Marr-Poggio  stereo  .algorithm.  This  al¬ 
gorithm  calculates  a  sparse  set  of  depth  values  from  the  disparities  of  all  the 


matched  zero  crossings  in  stereo  images.  Grimson  was  faced  with  the  problem 
of  interpolating  a  dense  surface  from  a  sparse  set  of  depth  measurements.  The 
problem  is  that  there  are  an  infinite  number  of  surfaces  which  pass  through  the 
measured  points.  Which  of  these  surfaces  is  most  likely  to  be  the  original  surface 
from  which  the  original  images  were  taken? 

Grimson  reasoned  that  there  should  be  no  wild  undulations  in  the  surface  in 
between  measured  depth  points.  If  wild  undulations  were  present,  then  it  would 
be  likely  that  more  points  would  be  marked  by  zero  crossings  and  produce  depth 
measurements  in  the  stereo  module.  He  called  this  the  “no  news  in  good  news” 
constraint.  Because  there  is  no  unique  solution,  surface  interpolation  is  an  ill- 
posed  problem  in  the  sense  of  Hadamard,  and  therefore  must  be  “regularized” 
[Poggio  and  Torre  1984'.  The  “no  news  is  good  news  constraint”  leads  to  the 
necessary  stabilizing  functional  which  ensures  the  uniqueness  of  the  solution. 

Using  this  smoothness  constraint  of  Grimson’s,  we  present  a  very  simple 
treatment  of  the  mathematics  for  the  two-dimensional  surface  interpolation  case. 
Surface  interpolation  requires  finding  a  surface  which  exactly  passes  through 
sparsely  measured  depth  points.  These  measurements  are  exact  and  no  other 
types  of  measurements  are  included  in  the  final  interpolated  surface.  A  possible 
set  of  such  measurements  are  depicted  in  Figure  2-2.  We  want  to  find  a  surface 
which  passes  exactly  through  these  points. 

Another  way  of  stating  the  “no  news  is  good  news  constraint"  is  that  we 
would  like  to  minimize  the  first  derivative  between  the  data  points.  Given  sparse 
depth  points  u  in  the  x-y  plane,  we  would  like  to  minimize  the  following  energy 
functional: 

£,=//(i,!  +  (5y)2‘'I<i!'  (21) 

The  energy  E\  can  be  thought  of  as  the  energy  of  a  rubber  membrane  'Courant 
and  Hilbert  1953.  p.  247j.  This  functional  should  only  be  minimized  in  the 


regions  between  the  depth  points.  The  exact  values  of  the  solution  surface  at 
the  depth  points  arc  fixed. 

We  now  use  the  Eulcr-Lagrangc  equations  [Con rant  and  Hilbert  1962|  to  find 
an  expression  for  the  stationary  point  where  the  minimum  energy  must  occur. 
An  application  of  the  Euler  Lagrange  equations  to  (2.1)  yields  the  following 
solution: 

Vau  =  0  (2.2) 

3  3 

in  between  the  constraint  points  (where  Va  =  (J^)  +  (^)  ).  This  is  equivalent 
to  fitting  a  flexible  membrane  between  the  constraint  points.  The  membrane 
solution  to  the  example  problem  is  shown  in  Figure  2-3. 

The  membrane  solution  can  also  be  modeled  with  a  mesh  of  resistors.  Here, 
the  voltages  are  proportional  to  depths  in  the  viewed  surface.  The  energy  term, 
equation  (2.1)  corresponds  to  the  energy  dissipated  in  all  of  the  resistors  in 
the  mesh.  Equation  (2.2)  is  a  restatement  of  KirchofT’s  current  law  that  all 
the  currents  flowing  into  a  node  must  sum  to  sero.  Constraints  are  added  by 
applying  voltage  sources  sparsely  throughout  the  mesh. 


Figure  2-3.  Membrane  Solution. 


There  arc  three  fundamental  problems  with  the  membrane  solution  of  sur¬ 
face  reconstruction: 

1.  As  seen  in  Figure  2-3,  the  membrane  does  not  cxtrapolaLc  beyond  the  con¬ 
straint  points.  The  surface  is  always  Rat  outside  the  given  points. 

2.  A  membrane  does  not  have  enough  rigidity  to  incorporate  orientation  con¬ 
straints. 

3.  It  has  been  shown  with  random  dot  stereograms  that  huinans  do  not  in¬ 
terpolate  surfaces  by  fitting  straight  lines  in  this  way  [Grimson  1981,  p. 
102].  We  tend  to  fit  smoother  surfaces  that  do  not  bend  sharply  at  the  data 
points  ,  for  example  see  Figure  2-4. 

A  more  realistic  interpolation  operator  must  be  derived.  A  better  operator 
results  from  an  energy  functional  which  minimises  the  sharp  bends  in  the  solu¬ 
tion,  minimising  the  change  in  the  first  derivative.  This  means  minimising  the 
second  derivative,  so 


Grimson  calls  this  energy  the  quadratic  variation  of  the  surface.  It  can  also  be 
modeled  as  the  energy  of  a  thin  plate  [Courant  and  Hilbert  1953].  Applying 


Figure  2-4.  Thin  plate  solution.  V 

Ruler’s  equations  to  minimise  Z?2»  we  get 

V2V2u  =-•  0  (2.4) 

which  is  equivalent  to  fitting  cubic  splines  in  each  interval.  Solution  or  this 
equation  yields  the  curve  surface  shown  in  Figure  2-4.  Notice  that  extrapolation 
is  done  around  the  exterior  and  the  surface  changes  smoothly. 

Crimson  developed  an  iterative  algorithm  based  upon  minimising  the  quadratic 
variation  lie  worked  on  the  surface  approximation  problem  since  he  did  not 
deal  with  surface  discontinuities  or  orientation  constraints.  Grimson’s  algorithm 
was  also  plagued  by  an  extremely  slow  convergence  rate. 


The  primary  contribution  of  Tcrsopoulos  [1984]  was  improved  computa¬ 
tional  efficiency  of  surface  reconstruction.  He  used  multigrid  methods  [Drandt 
1977]  to  speed  up  his  algorithm  by  orders  of  magnitude  over  the  single-grid 
approach.  Multigrid  speedup  ideas  will  be  discussed  in  detail  in  Chapter  4.  Ter- 
sopoulos  addressed  the  full  surface  reconstruction  problem  since  he  integrated 
both  depth  and  orientation  constraints  provided  by  various  visual  modalities.  In 
•addition,  he  studied  techniques  for  the  detection  and  the  explicit  representation 
of  surface  discontinuities.  As  mentioned  in  the  beginning  of  this  chapter,  the 


Figure  2-5.  21)  Itihnruionic  operator. 


location  of  discontinuities  is  assumed  to  be  input  to  the  surface  reconstruction 
module.  Some  the  .gills  about  how  we  might  combine  discontinuity  detection 
with  the  coupled  dcplh/slope  network  arc  discussed  iu  Chapter  6. 

Both  Crimson  and  Tcrzopoulos  used  the  standard  discrete  biharmonic  op¬ 
erator  for  their  surface  reconstruction  algorithms.  Crimson  used  Guile  difference 
methods  [1981,  p.  180)  to  arrive  at  the  biharmonic  mask.  Tersopoulos  used  more 
general  finite  element  techniques  but  chose  simple  uniform  squares  to  produce 
the  same  biharmonic  mask.  The  biharmonic  mask,  pictured  in  Figure  2-5,  can 
be  viewed  as  a  constraint  equation  relating  neighboring  depth  points.  The  mask 
can  equivalently  lie  expressed  as 

20u,7 

-  8(u,_ij  -f  +  u -f-Uij  +  i) 

f  2(«t  -  i.j- 1  F  Ui-ij  +  i  +  +  tti+ij  +  i)  (2.5) 

+  +  Utj- a  +  «.,/  ta)  =  0 

The  large  number  of  relations  of  this  type  form  an  enormous  system  of  linear 
equations  to  be  solved.  Many  solution  techniques  arc  available  for  solving  line<ur 

I 

systems  and  they  fall  logically  into  two  main  categories:  direct  and  iterative  J 

methods.  Direct  methods  take  a  Gnite  number  of  steps  to  come  up  with  an  I 

t 

answer.  Iterative  methods  start  with  some  initial  guess  of  the  solution  Mid  j 


gradually  converge.  Iterative  methods  arc  preferred  for  large  problems  such 
as  this  because  they  arc  faster,  cosier  to  program  and  take  less  storage  than 
the  corresponding  direct  methods.  Crimson  used  a  gradient  projection  iterative 
algorithm  while  Terzopoulos  chose  the  Gauss-Seidcl  iterative  method  on  each 
grid  in  a  hierarchy  of  grids  to  speed  convergence. 


The  Gauss-Seidcl  method  is  the  sequential  replacement  of  each  point  in  the 
mesh  by  some  function  of  other  points  in  the  mesh.  One  Gauss-Seidcl  realization 
of  (2.5)  is: 
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-2(V/,,.  i  d  ui  lj  +  l  +  ui i|+lj-  i  +  ttfnj  +  l) 

—  (u«  -~2j  +  Ui  12  J  +  Uf.y-2  +  ui.j+ 2)] 

Notice  that  each  node  replaces  itself  by  a  function  of  its  neighbors  some  of 
them  old  and  some  of  them  new.  For  this  reason,  the  Gauss-Seidcl  scheme  is 
inherently  sequential.  A  Jacobi  iterative  scheme  calls  for  the  replacement  of  all 
node  values  simultaneously,  and  so  can  be  parallel.  The  simple  extension  of  the 
Gauss-Seidel  relaxation  to  the  fully  parallel  Jacobi  scheme  results  in: 

‘  =  ^  W  +  <y-i  +  <y+i) 


-2(uf_ij_|  -f  u,fc_lj  +  J  -f  +  «,k+i,y+t)  (2.6) 

-Kfc~2j  +  2j  1-  uf./  2  +  «S,-+ 2)] 

Unfortunately,  the  standard  biharmonic  mask  is  unstable  when  used  in  a  par¬ 
allel  Jacobi  iteration.  See  Figure  2-6  for  an  example  of  this  instability.  The 
mathematical  reasons  for  the  instability  of  the  parallel  biharmonic  mask  are 
well  known.  See  [Varga  1962;  Young  1971 1  for  more  detail. 


flrandt  [1977]  describes  a  technique  called  the  weighting  of  corrections  to 
obtain  convergence  for  the  parallel  biharmonic  operator.  A  single  iteration  con¬ 
sists  of  two  steps:  a  biharmonic  application  and  a  local  weighting.  In  the  first 
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step  we  calculate  the  required  correction,  Sukj,  given  by  the  Jacobi  biharmonic 
replacement  (2.6),  where  uk  +  6k  =  u*+1.  That  is, 

6ukj  =  -uktj  +  ^  [8(uf_1(J  +  uf+li>  +  «?j_!  +  »,;,•+,) 

"h  **i— lj'+I  "h  tti+lj-l  u»+lj  +  l) 

-(«f-2,y  +  “*+2,y  +  «*;-2  +  tt?,y+2)] 

We  do  not  use  6v.ij  as  the  correction  term  since  the  method  does  not  converge. 


Instead  we  choose  some  linear  combination  of  the  local  6ut|,  as  the  correction 


term.  In  particular,  Brandt  suggests: 


*  =  *ij  +  1  j  +  +  6uiJ-l  +  6uiJ+ 1) 


For  the  optimal  smoothing  rate  Brandt  derived: 


u)a  =  1.552 


up  =  0.353 

Brandt  points  out  that  this  technique  is  very  sensitive  to  changes.  For  example, 
if  ua  were  changed  to  1.4  the  method  would  not  converge.  In  one  sense,  this 


weighting  scheme  is  equivalent  to  the  application  of  a  25-point  operator.  That 
is,  a  single  iteration  step  at  each  point  in  the  mesh  requires  the  values  of  24 
neighboring  points. 


There  is  no  natural  way  to  integrate  orientation  constraints  short  of  going 
back  to  the  original  energy  functional  and  adding  some  kind  of  penalty  term. 
In  fact  this  is  the  method  used  by  Terzopoulos  [1984].  At  points  where  the 
orientation  is  known,  he  creates  a  new  energy  term  which  is  the  square  of  the 


difference  between  the  orientation  measurement  and  the  explicitly  calculated 
slope.  Terzopoulos  adds  this  term  into  his  energy  functional  with  some  appro¬ 
priately  chosen  weighting  factor. 


Chapter  3  develops  the  coupled  depth/slope  model.  This  new  approach  to 
surface  reconstruction  leads  to  a  digital  algorithm  which  has  a  smaller  mask  size 
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than  the  weighted  technique  discussed  above.  Also,  orientation  and  discontinuity 
information  are  straightforwardly  integrated  into  the  surface  description. 
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CHAPTER  3 

THE  COUPLED  DEPTH/SLOPE  MODEL 


The  coupled  depth/slope  model  is  a  network  of  ideal  subtractor  elements 
connected  by  two  planes  of  resistor  meshes.  This  novel  model,  based  upon 
analog  network  considerations,  has  significant  advantages  over  previous  models 
of  surface  reconstruction. 

S.l.  ID  Form  of  the  Coupled  Depth/Slope  Model 


For  simplicity,  we  first  discuss  the  model  in  one  dimension.  The  proposed 
network  consists  of  depth  nodes  connected  to  slope  nodes  by  ideal  subtract 
constraint  boxes.  The  continuous  ID  network  is  shown  in  Figure  3-1.  A  resistor 
mesh  is  used  to  smooth  the  slope  points.  The  voltage  at  the  u  nodes  represent 
actual  depths  and  the  the  values  of  p  correspond  to  slope.  Note  that  u  and  p  are 
continuous  functions  of  z.  This  model  is  called  the  coupled  depth/slope  model 
because  of  the  coupling  between  the  depth  and  slope  representations  provided 
by  the  ideal  subtractor  elements.  The  subtractors  explicitly  calculate  a  slope 
representation  of  the  surface.  We  smooth  out  any  sudden  changes  in  the  slope 
with  a  resistor  mesh  and  pass  the  information  back  up  to  the  depth  plane.  Any 
depth  or  slope  node  can  be  made  into  a  constraint  by  fixing  an  ideal  voltage 
source  to  the  proper  place  in  the  network.  For  now,  assume  that  the  data  are 
exact,  forcing  the  reconstructed  depth  and  slope  values  to  pass  through  the  given 
constraint  values.  Noisy  and  conflicting  data  are  discussed  in  section  3.3. 


The  subtractor  device  is  shown  in  Figure  3-2.  If  nodes  A  and  B  are  set  with 
ideal  voltage  sources  then  node  C  will  be  set  to  A  -  B  by  the  device.  The  device 


Figure  3-1.  ID  network  for  coupled  depth /slope  surface  reconstruction 


is  unusual  in  that  all  of  its  terminals  can  act  as  inputs  or  outputs.  If  nodes  B 
and  C  arc  held  constant  with  voltage  sources,  then  the  A  terminal  is  fixed  to 
B  +  C.  If  A  and  C  arc  input,  then  B  becomes  A  —  C. 

Now  we  show  that  this  coupled  model  has  the  right  properties  for  surface 
reconstruction  as  dcGncd  in  Chapter  2.  The  model  must  Gnd  a  smooth  surface 
which  (its  the  given  depth  constraints.  In  addition,  orientation  constraints  and 
surface  discontinuities  must  be  incorporated  into  the  final  surface  description. 

First,  we  show  that  the  network  provides  for  biharmonic  smoothness,  used 
by  Grimson  and  Tcrzopoulos,  when  depth  constraints  only  arc  provided.  This  is 
shown  by  first  writing  the  power,  E,  dissipated  in  the  network.  The  minimization 
of  this  power  yields  the  biharmonic  equation.  The  power  dissipated  in  a  single 
resistor  equals  the  square  of  the  voltage  drop  through  the  resistor  times  the 
conductance  of  the  resistor  {E  =  gV7).  We  assume  that  the  vertical  resistors 
have  unity  conductance  and  the  ideal  sublractor  elements  consume  negligible 
power.  To  Gnd  the  total  power  dissipated  we  integrate  over  all  i  values. 

E  =  /(^-p)J  +  j(z)(^)3<<x  (3.1) 

The  first  and  second  terms  in  the  integral  represent  the  power  dissipated  in  the 
vertical  and  horizontal  resistors  of  the  network,  respectively. 


Maxwell’s  minimum  heat  theorem  states  that  the  distribution  of  currents 
and  voltages  in  a  circuit  is  such  that  the  total  power  dissipated  as  heat  is  mini¬ 
mized.  The  Euler- Lagrange  equation  are  used  to  find  the  minimum  energy  state 
reached  by  the  network.  The  appropriate  ID  Eulcr-Lagrange  equations  are: 

-£»  =  0  (3.2) 

—  —  Ep  —  0  (3-3) 

Assume  for  now  that  g(x)  =  1  for  .all  x.  Applying  these  Eulcr-Lagrange  equations 
to  find  the  stationary  points  of  the  energy  functional,  E,  results  in: 

s*  =  t  <3<> 

£  =  «■-=  <3-5> 

These  are  a  pair  of  coupled  second  order  equations  relating  u  and  p.  These 
equations  are  shown  to  reduce  to  the  bihannonic  equation  by  first  taking  of 
(3.4)  and  of  (3.5)  giving: 

£  -  £  <■  *> 

S-i-a  «") 


*  «l*  «  .  r,  »  ,r,  ,  *  .  •  .»  .*  «  .  r  ,  t  ,  **■*•  ir  W,  n*  r ,  v.  r, 


24 


or 


dx 


dx2 


d*u  dp  d2u 

(3.8) 

dx 4  dx  dx 2 

=  0  from  (3.4),  so 

d*u 

(3.9) 

d^=° 

This  is  the  one-dimensional  form  of  the  biharmonic  equation. 

3.2.  Similarities  to  the  Human  Visual  System 

There  is  a  lot  of  redundancy  in  the  2^-D  sketch  described  here  because  both 
depth  and  slopes  are  explicitly  represented.  Only  one  of  these  measures  need 
by  stored,  and  the  other  can  be  recovered  with  a  simple  derivative  or  integral. 
Marr  '1982,  p.  279]  draws  on  psychophysical  evidence  to  argue  that  the  human 
biological  2|-D  sketch  makes  explicit  representations  of  orientation  as  well  as 
depth. 

3.3.  Extending  the  Model  to  Handle  Noisy  Constraints 

Thus  far  we  have  assumed  that  the  given  depth  and  orientation  constraints 
are  exact.  We  force  the  depth  and  slope  planes  to  pass  exactly  through  the 
constraint  points.  In  reality,  there  are  additive  noise  components  for  each  con¬ 
straint  measurement.  It  would  be  unwise  to  force  the  reconstructed  surface  to 
pass  through  noisy,  inaccurate  points.  This  problem  is  solved  by  adding  another 
term  to  the  energy  functional.  We  rewrite  the  energy  in  (3.1)  as: 

E  =  J (u *  -  p)2  +  (Pl)2  +  a(z)(u  -  u)2  +  0{x){p  -  p)2  dx  (3.10) 

The  u  and  p  are  the  values  of  the  depth  and  slope  measurements  respectively. 
The  variables  a(x)  and  3(x)  represent  the  confidence  of  the  measured  values  in 
depth  and  slope.  For  example,  if  we  knew  that  the  measurement  u(xo)  was 
inaccurate  we  could  set  the  corresponding  a(io)  to  zero. 
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(3.11) 

(3.12) 


p(x)  uricnliition  constraint 

Figure  3-3.  Handling  noisy  constraints. 

Similar  energy  terms  were  used  by  Tcrzopotilos,  Crimson  and  Marroquin. 
According  to  Marroquin  [1985], 

“W  =  2o2(u) 

= s4)  <3,2> 

where  o(u)  represents  the  standard  deviation  of  zero  mean,  white,  additive  Gaus¬ 
sian  noise  of  the  measurement  u.  Therefore  a(x)  and  (3[x)  represent  a  tradeoff 
between  surface  smoothness  and  how  close  the  surface  Gts  to  the  depth  and 
orientation  constraints. 

The  added  noise  term  can  be  modeled  as  a  conductance  times  the  square  of 
a  voltage  and  so  is  easily  added  to  our  electrical  model.  We  now  set  constraints 
with  a  voltage  source  through  a  series  resistor  to  the  constraint  point  in  the 
network.  These  connections  arc  shown  in  Figure  3-3.  The  voltage  source  has  a 
value  of  the  depth  measurement  u.  This  voltage  is  applied  through  a  resistor 
with  conductance  ar(x),  where  a(x)  is  related  to  the  amount  of  noise  in  u  by 
equation  (3.11).  If  there  is  no  depth  constraint  at  a  particular  location,  say  xo, 
we  choose  an  arbitrary  constraint  value  u(xq)  for  the  voltage  source  with  an 
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open  circuit,  or(zo)  =  0.  This  is  equivalent  to  saying  that  there  is  infinite  noise 
in  u(xo)  and  that  we  have  no  confidence  in  that  surface  measurement. 


Minimizing  (3.10)  gives: 


<f2tt  dp  I  .X 
d7p  du 


(3.13) 

(3.14) 


The  solution  of  these  equations  gives  a  smooth  surface  which  perforins  a 
least  squares  averaging  reduction  of  noise.  A  very  different  method  for  handling 
noise  is  proposed  in  Chapter  4.  In  a  digital  algorithm,  we  may  want  to  assume 
noiseless  measurements  and  solve  the  problem  exactly.  Then  we  can  perform 
some  amount  of  smoothing  dependent  on  how  much  noise  is  present.  This  post 
smoothing  step  is  very  fast  since  we  are  dealing  with  dense  surface  data. 

3.4.  Generalized  Smoothness 

How  much  smoothness  should  we  require  of  a  surface?  Some  surfaces  are 
very  jagged  and  the  biharmonic  reconstructs  a  surface  that  is  much  too  smooth. 
On  the  other  hand,  some  surfaces  are  even  more  smooth  than  the  biharmonic; 
they  may  have  smooth  second  or  third  derivatives.  The  biharmonic  is  typically 
used  for  surface  interpolation  because  it  seems  to  do  reasonably  well  over  all 
possible  surfaces. 

Terzopoulos  [1986]  proposes  Tikhonov  stabilizers  [Tikhonov  1977]  to  gen¬ 
eralize  smoothness  to  any  order  of  derivative.  For  example,  the  ID  Tikhonov 


stabilizer  appropriate  for  our  problem  is: 


E=  £  /  Srn(*)(~^)3dz 


(3.15) 


where  the  <7m(x)  are  continuous  weighting  functions  and  is  the  m,fl  deriva¬ 

tive  of  u.  The  generalized  coupled  depth  /slope  network  directly  implements 
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the  control  left-continuity  stabilizers  of  Terzopoulos.  The  generalized  network  is 
shown  in  Figure  3-4.  Here  uu(i)  represents  the  depth  values,  while  «i(i)  repre¬ 
sents  the  first  derivative  of  uo(x).  Naturally,  um(z)  represents  the  mth  derivative 
function  of  uo(x),  and  j„,(z)  denotes  the  conductance  function  of  level  m  of  the 


network. 


As  before,  we  can  write  the  energy  dissipated  in  all  the  resistors  in  the 


network  as: 


E  =  £  / +  l^*.W  -  .,(*))’]<<*  (3.16) 


The  first  term  in  the  integral  represents  the  power  dissipated  in  the  horizontal 
resistors  on  the  m<fc  level,  .and  the  second  term  is  the  power  lost  in  the  vertical 


resistors. 


This  very  general  network  can  incorporate  constraints  of  .any  order  of  deriva¬ 
tive  into  the  surface  solution.  This  is  done  simply  by  setting  the  proper  mesh 
points  with  a  voltage  source.  Also,  discontinuities  in  the  m<fc  order  derivative  at 
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Figure  3-5.  A  depth  discontinuity. 


i  —  i0  can  be  incorporated  by  setting  ffm(* o)  =  0.  For  example,  a  depth  discon¬ 
tinuity  can  be  implemented  by  such  a  circuit  as  shown  in  Figure  3-5.  Here,  we 
have  broken  the  resistors  that  smooth  the  first  (and  higher)  derivatives.  Depth 
discontinuities  occur  in  between  points  in  the  u  mesh. 

Terzopoulos  handles  orientation  discontinuities  by  switching  to  a  membrane 
model  instead  of  the  thin  plate  at  the  location  of  the  discontinuity.  This  allows 
the  reconstructed  surface  to  crease  sharply  but  still  remain  continuous.  In  the 
coupled  dcpth/slopc  model,  membrane  smoothness  is  provided  by  the  top  re¬ 
sistor  mesh,  which  smoothes  the  depth  values.  The  resistors  on  all  other  levels 
are  broken.  Orientation  discontinuities  occur  exactly  on  u  mesh  locations.  An 
example  of  an  orientation  discontinuity  is  shown  in  Figure  3-6.  The  thin  plate 
mode!  is  not  rigid  enough  to  allow  other  even  higher  order  discontinuities,  but 
the  coupled  dcpth/slope  network  can  handle  discontinuities  of  any  order. 


orientation  discontinuity 


Figure  3-0.  An  orientntion  discontinuity. 

For  some  surfaces,  the  thin  plate  solution  is  too  smooth.  However,  we  cannot 
use  the  membrane  solution  because  of  its  tendency  to  crease  and  its  poor  ex¬ 
trapolation  properties.  Splines  under  tension  were  invented  to  combat  a  similar 
problem  in  graphics  [Pilcher  1974],  Splines  under  tension  provide  a  mechanism 
for  combining  varying  amounts  of  smoothness  from  the  thin  plate  and  mem¬ 
brane  smoothness.  They  are  also  proposed  as  a  surface  interpolation  technique 
to  reconstruct  surfaces  with  intermediate  smoothness  between  thin  plates  and 
membranes  [Terzopoulos  1984].  The  coupled  dcpth/slopc  network  generalizes 
splines  under  tension  since  the  network  can  integrate  .arbitrary  combinations  of 
membrane,  thin  plate  and  any  higher  order  smoothness. 

The  general  coupled  dcpth/slope  network  generalizes  the  surface  reconstruc¬ 
tion  algorithm  of  Terzopoulos  [1984].  Information  at  arbitrary  orders  of  smooth¬ 
ness  can  be  usefully  combined  into  the  surface  description. 


Figure  3-7.  2D  extension  -  the  p  plane. 


3.5.  2D  Extension  of  the  Coupled  Dcpth/Slope  Method 

The  2D  network  for  the  coupled  dcpth/slope  method  is  shown  in  Figures 
3-7  and  3-8.  For  simplicity,  we  revert  back  to  biharmonic  smoothness  only,  but 
the  generalized  smoothness  ideas  developed  in  the  last  section  apply  equally  well 
in  two  dimensions.  The  2D  network  is  more  complicated  because  the  slopes  in 
both  the  x  and  the  y  directions  must  be  included.  Therefore  we  use  two  resistor 
meshes  to  smooth  the  p  and  q  slope  planes,  where  p  is  the  slope  of  u  in  the  x 
direction  and  q  is  the  slope  in  the  y  direction. 

Again,  the  power  dissipated  in  the  total  network  is  found  to  be: 

E  =  j  J  [(«*  -  P)2  +  K  ~  9)2  +  P*  +  pJ  +  9*  +  ?j]  dxdV  (317) 

The  extension  for  noisy  constraints  described  in  section  3.3  for  the  ID  case  is 
also  applicable  here.  The  appropriate  2D  Euler-Lagrangc  equations  for  (3.17) 
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V2u  =  px+qy 

(3.18) 

V2p  =  p-uz 

(3.19) 

V2q  =  q  -  ut 

(3.20) 

w^re  V*  =  (*)*  +  (*)*- 

Taking  V2  of  (3.18)  and  of  (3.19)  and  ^  of  (3.20)  gives: 

V2V2«  =  V2pt  + 

V2p,  =  p»  -  »«« 

=  9*  ~  •**» 

Substituting  (3.22)  and  (3.23)  into  (3.21)  yields: 

V2V2u  =  pz  +  qy  -  uXi  -  uvy  =  Px  f  9V  -  V2u 
We  know  from  equation  (3.18)  that  p»  +  9V  -  V2ti  =  0,  so 

V2V2u  =  0 

which  is  the  two-dimensional  biharmonk  equation. 


(3.21) 

(3.22) 

(3.23) 


(3.24) 


From  the  coupled  differential  equations,  we  cannot  show  that  p  and  q  must 
converge  to  uz  and  u„.  A  simple  combination  of  the  original  three  coupled 
equations  gives: 


V2V2u 


d  d 


We  can  only  show  that  uz  =  p  and  uv  ~  q  if  some  boundary  conditions  are 
assumed  or  else  we  define  pv  —  qz. 


3.6.  Other  Computational  Models 


Computational  models  which  solve  the  biharmonic  equation  must  neces¬ 
sarily  rontain  some  active  devices.  It  can  be  shown  that  the  biharmonic  mask 
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Figure  3-9.  Arbitrary  configuration  of  resistors  and  sources. 

cannot  be  synthesized  out  of  purely  positive  valued  resistive  components.  The 
extrapolation  mechanism  required  for  biharmonic  smoothness  can  boost  the  volt¬ 
age  on  the  exteriors  of  the  array  to  higher  and  lower  voltages  than  that  which 
appear  interior  to  the  array.  The  resistor  network  alone  cannot  boost  a  voltage 
at  a  node  higher  than  the  voltage  of  any  of  the  sources  in  the  system.  This  is 
seen  in  Figure  3-9.  The  voltage  at  node  A,  for  example,  cannot  be  higher  than 
the  voltage  at  its  neighboring  nodes.  If  this  were  the  case,  than  all  currents  will 
flow  out  of  node  A,  violating  KirchofT’s  Current  Law  (KCL). 

Nonetheless,  there  are  applications  for  the  simple  membrane  solution  given 
by  a  mesh  of  equal  valued  resistors.  For  example,  controlled  industrial  environ¬ 
ments  consisting  of  flat,  zero-slope  planes,  problems  with  undetectable  disconti¬ 
nuities,  and  scenes  with  no  discontinuities  arc  nil  good  applications  for  membrane 
smoothness.  However,  we  should  not  attempt  to  tackle  the  full  vision  problem 
without  a  more  rigid  smoothing  technique.  In  general,  the  membrane  solution 
gives  a  hist,  first  approximation  to  the  surface. 

Another  way  to  solve  the  biharmonic  equation  is  to  split  it  into  two  Poisson 
equations  and  solve  them  sequentially.  It  will  be  shown  that  this  procedure  only 
results  in  an  approximate  solution  to  our  surface  interpolation  problem. 


For  example,  the  bihannonic  equation  is  given  by: 


V4u  =  /  h-  V2V2u  =  / 

The  two  poisson  equations  are  as  follows: 

(3.25) 

V2v  =  / 

(3.26) 

V2u  =  v 

(3.27) 

We  solve  for  v  in  (3.26)  and  then  use  the  result  to  solve  for  u  in  (3.27). 


An  actual  implementation  of  this  method  shows  that  the  solution  of  the 
first  equation  is  slow  compared  to  the  solution  of  the  second.  This  is  because 
we  are  solving  a  very  sparse  system  in  the  first  equation.  The  solution  of  the 
second  equation  merely  results  in  local  perturbations  of  the  dense  values  of  v. 
Since  the  values  of  v  arc  equally  spaced,  we  could  solve  the  second  equation 
with  a  simple  gaussian  convolution  step.  The  convolution  method  has  been 
implemented  since  it  has  an  additional  advantage:  any  amount  of  smoothness 
can  be  integrated  into  the  surface  by  varying  the  siac  of  the  Gaussian  operator. 
A  very  narrow  Gaussian  will  have  little  effect  on  the  surface  and  so  return 
membrane  smoothness.  A  wider  Gaussian  will  return  thin  plate  and  smoother 
surfaces.  Using  either  method,  the  time  it  takes  to  solve  the  first  equation  is 
still  much  greater  than  the  time  to  solve  the  second. 

Thus,  the  dual  Poisson  method  is  another  way  to  generate  a  quick,  crude 
solution  to  the  surface  interpolation  problem.  This  method  improves  upon  the 
membrane  solution  by  forcing  the  surface  to  be  smooth  without  too  much  extra 
computation.  The  main  difficulty  with  the  dual-Poisson  approach  is  the  bound¬ 
ary  conditions.  Either  the  slope  or  the  depth  values  must  be  specified  around 
the  boundary. 

3.7.  Summary 

In  summary,  we  have  developed  the  coupled  depth/slope  model  for  surface 
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reconstruction.  This  is  a  computational  model  which  has  a  number  of  features: 

1.  If  depths  only  are  given,  the  model  reconstructs  a  surface  with  the  same 
biharmonic  smoothness  used  by  Terzopoulos  and  Grimson.  In  addition,  in 
the  process  of  computing  the  depths,  a  dense  slope  representation  is  also 
calculated.  For  this  reason,  orientation  constraints  can  be  incorporated  into 
the  surface  description  with  no  extra  work. 

2.  Like  Terzopoulos’  model,  discontinuities  in  depth  and  orientation  can  be 
used  to  compute  the  reconstructed  surface.  However,  the  coupled  depth/slope 
model  generalizes  to  handle  constraints  of  any  order  derivative.  In  addition, 
the  mode]  can  handle  discontinuities  at  any  order  of  derivative  of  the  sur¬ 
face. 

3.  The  model  seems  to  be  biologically  feasible.  It  fits  the  four  criteria  given  in 
Chapter  1:  parallelism,  uniformity,  locality  and  fault  tolerance.  The  model’s 
biharmonic  smoothness  properties  are  observed  in  psychophysical  tests  of 
humans.  Also,  Marr  argues  that  the  human  2^-D  sketch  redundantly  repre¬ 
sents  both  slopes  and  depths.  This  compares  favorably  with  the  dense  depth 
and  slope  representations  generated  by  the  coupled  depth/slope  network. 

4.  Finally,  both  digital  and  analog  implementations  of  the  coupled  depth /slope 
method  are  possible. 
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CHAPTER  4 

DIGITAL  IMPLEMENTATION 


This  chapter  develops  a  digital  algorithm  based  upon  the  coupled  depth/slope 
model.  Again,  we  start  with  the  ID  case  before  we  progress  to  the  2D  ver¬ 
sion.  We  apply  standard  finite  difference  methods  to  the  continuous  differential 
equations  derived  in  the  last  chapter  to  get  a  large  system  of  linear  equations, 
which  we  solve  with  a  parallel  iterative  algorithm.  It  is  shown  that  this  itera¬ 
tive  algorithm  is  very  well  suited  to  a  mesh-of-processors  computer  such  as  the 
Connection  Machine. 


4.1.  The  Discrete  ID  Case 


The  coupled  differential  equations  for  u  and  p  derived  in  Chapter  3  were: 

d2  u  _  dp 
dx 2  dx 

d2p  du 


dx 2  ^  dx 


(3.4) 

(3.5) 


Applying  simple  finite  difference  approximations  for  first  and  second  derivatives 
[Smith  1978,  p.  6],  as  illustrated  in  Figure  4-1  yields: 


«i_i  -  2u,  +  u.+i  =  pi  -  pi- 1 

Pi- 1  -  2 pi  +  p,+ 1  =  pi  -  (u,+  i  -  u<) 


Rearranging  to  solve  for  u,  and  p,  gives: 


tit  =  -  [(**• — 1  +  Pi- 1)  +  (ti»+l  -  P«)] 

(4.1) 

Pt  =  -  l(u,+i  -  tit)  +  Pt-I  +  Pt  +  l] 

(4.2) 

A  full  Jacobi  iteration  method  (simultaneous  displacement)  for  these  equa¬ 
tions  fails  to  converge  to  a  solution.  A  proof  of  this  divergence  for  the  2D  case  will 
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Figure  4-1.  Discrete  coupled  dcpth/slope  model. 


be  shown  in  section  4.2.  However,  a  simple  alternating  process  of  first  calculat¬ 
ing  the  new  depth  values  and  then  the  new  slope  values  guarantees  convergence. 


The  iteration  method  is  therefore: 


«‘+*  =  *•[«-.  +»£.)  +  (■*?+. -rf)] 
=  5[<“.‘+'  -  “?+1)  +  +  pf+.l 


where  the  superscript  k  denote  the  iteration  step.  A  proof  of  convergence  for 


the  2D  the  alternating  depth/slopc  algorithm  is  given  in  section  4.2. 


4.2.  The  2D  Discrete  Case 


We  would  now  like  to  extend  our  analysis  to  two-dimensions.  We  know 
from  Chapter  3, 

V2u  =  p, +  7,  (3.16) 


V2p  =  p  -  ii« 


(3.17) 


v’fl  =  q-u. 


(3.18) 


Applying  2D  finite  difference  approximations  to  the  coupled  differential  equa¬ 
tions  and  rearranging  terms  results  in  the  following  discrete  equations: 
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u«,j  =  +  Pi-i, j)  +  +  9ij- 1)  +  («.j  +  i  -  9«,j)  +  (u.+»j  -  P.,j)] 

Ptj  =  gI(U«+l.J  “  U*j)  +  Pi+lJ  +  Pi-1,}  +  Pij+1  +  PtJ-li 

9tj  =  g [(«,,>+ 1  -  «ij)  +  9»+i,j  +  9« — 1  j  +  9«,y+i  +  9.J-1:  (4-5) 

These  finite  difference  equations  make  up  a  large  but  sparse  system  of  equations 
that  must  be  solved.  Because  of  the  size  of  the  system,  iterative  methods  are 
generally  used.  Standard  iterative  techniques  will  be  used,  but  other  alternatives 
are  discussed  later  in  the  chapter. 
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The  discretized  equations  for  tt,  p  and  q  derived  above  map  directly  into 
a  Gauss-Seidel  scheme  for  a  conventional  sequential  computer.  This  direction, 
though  important,  is  not  developed.  We  will  concentrate  on  Jacobi  iterative 
schemes  which  are  suitable  for  fine  grain,  parallel  computers.  The  architectures 
of  primary  concern  are  locally  connected  meshes  of  very  simple  processing  ele¬ 
ments  such  as  the  Connection  Machine  [Hillis  1985*. 

These  architectures  have  a  number  of  features  in  common.  They  all  consist 
of  a  2D  rectangular  mesh  of  tens  of  thousands  of  extremely  simple  processing 
elements.  Inter-processor  communication  is  performed  through  nearest  neigh¬ 
bor  connections  only.  Machines  such  as  the  Connection  Machine  have  special 
hardware  for  arbitrary  global  communication  but  this  communication  is  much 
slower  than  local  communication  and  is  not  appropriate  for  any  of  the  algo¬ 
rithms  developed  in  this  report.  The  processing  elements  typically  contain  only 
a  one-bit  ALU.  This  necessarily  requires  that  all  arithmetic  be  performed  bit 
serially.  Since  many  processors  are  packaged  per  chip,  each  processing  element 
is  allotted  only  a  limited  amount  of  on  chip  memory. 
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The  algorithms  developed  in  this  chapter  have  been  devised  to  match  the 
architecture  constraints  of  machines  such  as  the  Connection  Machine,  so  that 
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surfaces  can  be  calculated  as  fast  as  possible.  Thus  the  algorithms  exhibit  local 
operations  which  can  be  performed  in  parallel  over  all  points  in  the  mesh,  can 
be  calculated  easily  with  fixed  point  arithmetic,  and  require  very  little  memory. 


A  first  approach  at  developing  a  Jacobi  iterative  scheme  can  be  constructed 
from  the  coupled  difference  equations  (4.5): 


+  Pi- lj)  +  +  *ij- 1)  +  (“.J+i  " 

PiJ1  ~  7  KU‘fc+l,J  —  U£>)  +  P<+ 1 J  +  Pi~ 1 J  P‘J  +  1  +  PiJ-l] 


1  =  -  u*j)  +  9i+ u  +  9,-- ij  +  9.-.J  +  1  + 


Like  its  ID  counterpart,  this  method  does  not  converge  to  a  solution.  This 
should  not  be  surprising  since  it  has  been  shown  that  this  method  solves  the 
biharmonic  equation  for  depths  and  that  the  Jacobi  biharmonic  iteration  does 
not  converge.  Local  Fourier  analysis  is  now  be  used  to  show  rigorously  that  the 


above  iteration  method  does  not  converge. 


First,  we  define  the  error  values  Au*,  Apfc,  and  Aqk  as  the  differences 
between  the  values  of  u,  p,  and  q  at  a  given  iteration  step  k,  and  the  final  values 


of  the  true  solution. 


£,  =  UiJ  ~  ukj 


APi,j  =  PiJ  ~  Px,} 


A<lij  =  ?!,/ 


It  is  hoped  that: 


lim  Auf,  =  0 


lim  A  pi ,  =  0 

k  ■—  oo  1,3 


lim  A q*j  =  0 

it  — oo 


Similarly  we  define  the  errors  between  the  values  at  step  k  +  1  And  the  final 


solution  as: 
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These  equations  can  be  expanded  by  taking  the  difference  of  equations  (4.5) 
from  (4.6): 

Au?r  =  ^[(Au*-U  +  AP?-ij)  +  1  +  1) 

+  (Att*j  +  i  ~  +  (A«i+ij  -  AP*j)] 

&Pij 1  =  “  Au*j) 

+  AP?+i,i  +  AP*-ij  +  AP?J  +  i  +  APfj-i] 

=  ^(A«*y+1  -  Au?u) 

+  A?.k+i,>  +  A?*_ i,y  *  Arfjvi  +  A*. j  - »J  (4-7) 

Now,  let’s  assume  that  the  error  occurs  with  some  Fourier  spatial  frequency 
(<vj,u>2)  over  the  surface.  We  assume  a  different  Fourier  error  for  both  the  depth 


and  slope  values. 


Au*tJ  =  Uke'liw'^jUi) 

Apfj  =.  Pkei(,u>  +)“*) 
A  qfj  =  Qke'('“‘ 


Ai»fc+1  =  Uk+ 1  «•(“*'« +>u,3) 

Ap*+l  =  Pfc+V  (iwi  +;u»3) 

+  1  =  Qfc+le*(tui+;w2) 

where  (/,  P  and  Q  are  the  magnitudes  of  the  associated  Fourier  error.  Sub¬ 
stituting  these  expressions  into  the  equation  (4.7),  dividing  by  ex^lu>,+:'u)3^  and 
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rearranging  terms  gives: 


4  0  0 


where 


0  5  0 
0  0  5 


y  y  \h 
y  a  o 


/  =  coa(uq)  +  cos(cv2) 
g  =  eiw ‘  -  1 


h  =  e* 


The  amplification  matrix  A(u  1,0/2)  shows  how  much  the  magnitude  of  the 
errors  decreases  for  each  iteration.  A(u>i,u/2)  is  defined  as: 

fu\k+x  fu\k 

P  =  A(u  1,0/2)  P 

<Qj  W> 

A  simple  calculation  for  cjj  =  wj  —  ir  yields  eigenvalues  greater  than  one 
in  magnitude.  This  means  that  certain  frequency  components  of  the  error  are 
amplified  during  each  iteration,  preventing  the  calculation  from  converging. 

However,  we  can  add  more  stability  to  the  iteration  by  requiring  that  we 
alternate  between  depth  and  slope  values.  That  is,  first  we  calculate  all  of 
the  depths  in  parallel  and  then  we  calculate  all  of  the  slopes  in  parallel.  This 
iteration  scheme  looks  like  the  following; 

=  ^  [(«*-!./  +P?-w)  +  Ki-i  +  1)  +  («*y+i  -  <7  ij)  +  (uf+ij  ~  P?,y)] 

P*y  1  =  ^[(u*+ij  -  «.*y  *)  +  +  Pf-J.j  +P?j+x  t-Po-il 

Qtj  ~  5IK.J  +  1  ~  u.,y  )  +  9,+ij  +  <7,-i,y  +  <7,,y+i  + 

Notice  that  the  new  ufc+1  values  arc  computed  using  old  values  of  depths  and 
slope  (ufe,  pk ,  and  qk).  However,  the  new  pk+I  and  q*+1  values  need  the  newly 
computed  values  of  the  depth  (ufc  fl).  There  is  an  inherent  sequential  nature  to 
the  algorithm  since  the  depths  must  be  computed  first  and  then  the  slopes. 


Going  through  the  same  type  of  error  analysis  as  above  yields  the  Fourier 
error  equation  of: 

/  4  0  0\  /U\k+1  (\f  \g  \h\  / U\k 

-g  5  0  P  =  0  \f  0  P 

\-h  0  5  \QJ  VO  0  If)  [qJ 

with  /,  g  and  h  defined  as  above.  This  results  in  an  amplification  matrix 

A(u>uLJ2)  of 

(  }/  i»  i*  \ 
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It  has  been  shown  numerically  with  the  EISPACK  routines  [Goos  and  Hartmanis 
1976]  that  the  maximum  eigenvalue  of  the  amplification  matrix  for  this  case  has 
a  value  of  .953  at  Ui  =  u>2  =  2.79.  This  shows  that  the  method  converges  but 
has  a  very  high  smoothing  rate  and  so  would  not  be  appropriate  for  a  standard 
multigrid  implementation. 

This  algorithm  fits  nicely  on  a  computer  such  as  the  Connection  Machine. 
Assume  for  now,  that  there  are  enough  processors  to  fit  each  utJ  value  in  each 
processor.  The  corresponding  ptJ  and  qti]  values  are  also  stored  in  the  same 
processor  that  contains  ut(J.  This  mapping  scheme  is  shown  in  Figure  4-3. 

There  are  two  key  advantages  of  this  scheme.  First,  we  pay  no  penalty  for 
sequentially  computing  all  of  the  depths  (in  parallel  of  course)  and  then  com¬ 
puting  all  of  the  slopes  (again  in  parallel).  We  couldn’t  compute  the  depths  and 
slopes  simultaneously  anyway  since  they  are  calculated  using  the  same  proces- 


Secondly,  the  iterative  scheme  requires  that  processors  communicate  only 
with  their  nearest  four  neighbors.  This  is  in  direct  contrast  to  a  13-point  bi¬ 
harmonic  mask  implementation,  where  12  neighboring  values  must  be  collected. 
To  make  the  13-point  biharmonic  operator  converge,  an  extra  local  weighting 
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Figure  4-3.  Mapping  onto  the  Connection  Machine 


step  must  be  added  which  turns  each  iteration  into  the  equivalent  of  a  25-point 
operator  application. 


Another  advantage  of  this  scheme  is  that  it  requires  no  floating  point  op¬ 
erations.  Floating  point  operations  take  thousands  of  cycles  on  processors  with 
single  bit  ALUs  typically  found  in  these  machines.  These  serial  ALUs  actually 
give  these  machines  an  advantage  over  conventional  16-bit  or  32-bit  computers. 
The  inputs  to  the  surface  reconstruction  module  are  fairly  inaccurate  and  noisy. 
The  most  accuracy  we  should  expect  is  about  8  bits.  It  is  advantageous  to  only 
have  to  do  8-10  bit  word  arithmetic  bit-scrially.  A  conventional  32-bit  computer 
gains  nothing  from  this  reduced  word  size. 


4.3.  Examples  of  Surface  Reconstruction 

The  coupled  dcpth/slopc  algorithm  has  been  implemented  and  tested  on 
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a  number  of  different  synthetic  surfaces.  Figure  4-4a  shows  some  sparse  data 
measurements  that  have  been  sampled  from  a  synthetic  surface.  The  surface  re¬ 
constructed  by  the  coupled  depth/slope  network  is  shown  in  Figure  4-4b.  Figure 
4-4c  shows  the  result  given  the  location  of  the  circular  discontinuity. 

The  next  example,  shown  in  4-5  shows  the  reconstruction  of  a  surface  given 
only  one  u  constraint  and  100%  of  the  p  constraints.  In  effect,  the  coupled  depth 
slope  algorithm  integrates  the  slope  values  and  uses  the  single  depth  constraint 
to  find  the  constant  of  integration. 

Ikeuchi  discusses  a  similar  type  problem  in  [Ikeuchi  1983].  Here,  some  shape 
from  shading  module  has  produced  the  p  and  q  slopes  everywhere  in  the  image. 
Ikeuchi  employs  the  following  least  squares  energy  functional  shown  to  recover 
the  surface: 

£  =  /  J (“I  ~P)7  +  («y  -  qfdxdy  (4.7) 

Minimizing  this  energy  gives, 

V2u  =  pz  +  qy  (4.8) 

This  leads  to  a  simple  iterative  scheme  to  recover  depth  from  the  slope  mea¬ 
surements.  It  should  be  noted  that  equation  (4.8)  derived  by  Ikeuchi  is  identical 
to  the  first  of  the  three  equations  (3.18)  describing  the  coupled  depth/slope 
method.  The  coupled  depth/slope  method  is  therefore  a  generalization  of  the 
depth  from  shape  algorithm  developed  by  Ikeuchi. 

How  many  constraints  must  be  given  in  order  for  the  coupled  depth/slope 
problem  to  be  well  formulated?  Terzopoulos  [1984,  p.  73]  gives  the  minimal  set 
of  input  constraints  for  the  surface  reconstruction  problem  to  be  well-posed  as: 

1.  Three  noncolinear  depth  constraints, 

2.  Two  depth  constraints  as  well  as  a  single  (nonaligned)  p  or  q  constraint, 

3.  A  single  depth  constraint  as  well  as  a  single  p  and  a  single  q  constraint. 
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4.4.  Extension  for  Handling  Noisy  Constraints 

We  now  propose  a  different  way  for  handling  noise  than  that  discussed  in 
Chapter  3.  In  a  digital  algorithm  on  a  Connection  Machine,  we  assume  noiseless 
measurements  and  solve  the  problem  exactly.  This  keeps  the  critical  loop  which 
determines  the  speed  of  the  algorithm  as  fast  and  simple  as  possible.  After  an 
exact  solution  is  reached,  we  do  some  amount  of  smoothing  dependent  on  how 
much  noise  is  present  in  the  constraint  measurements.  This  post  smoothing  step 
is  very  fast  since  it  deals  with  dense  data. 

An  example  of  this  post  smoothing  approach  for  handling  noisy  constraints 
is  shown  in  Figure  4-6.  Initially,  the  exact  solution  is  very  bumpy.  We  perform 
a  number  of  smoothing  steps  proportional  to  how  much  noise  we  believe  is  in 
the  data.  These  smoothing  steps  reduce  the  bumps  in  the  surface.  Different 
confidence  values  are  handled  by  allowing  the  high  noise  points  to  change  more 
than  the  low  noise  ones.  This  technique  is  similar  but  not  equivalent  to  the  least 
squares  averaging  approach  used  by  Crimson  and  Terzopoulos. 

4.5.  Fault  Tolerance 


m 


Fault  tolerance  is  a  significant  problem  that  has  to  be  dealt  with  in  massively 
parallel  architectures.  Suppose  we  fabricate  a  32  x  32  array  of  processing  ele¬ 
ments  on  a  wafer  and  that  10^  of  them  are  broken.  Nevertheless,  we  would  like 
the  algorithm’s  performance  to  degrade  gracefully  as  more  processors  become 
faulty.  A  diagnostic  routine  can  be  run  in  which  processors  check  their  neighbors 
to  see  if  they  are  functioning  properly.  If  a  neighbor  is  faulty,  a  discontinuity  is 
marked  in  that  direction. 


Notice  here  that  it  doesn’t  matter  if  a  neighboring  processor  is  faulty  or  if 
the  wire  connecting  them  is  broken.  For  example,  in  Figure  4-7,  processor  A 
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a)  exact  solution 


b)  after  1  smoothing  step 


c)  after  20  smoothing  steps 


Figure  4-0.  Handling  noise. 


marks  B  jus  being  broken  since  A  gets  no  response  over  their  connecting  wire.  B’s 


other  neighbors  however  think  that  B  is  in  fine  shape  and  can  do  a  good  job.  Of 
course,  B  thinks  that  it  is  A  that  is  broken  since  A  never  responds  to  his  queries. 
An  example  of  a  solution  of  a  problem  given  10%  random  broken  processors  out 
of  an  array  of  32  by  32  is  given  in  Figure  4-8.  Notice  that  the  overall  solution 
stays  very  much  correct,  we  only  lose  a  small  amount  of  resolution  at  the  location 
of  bad  processor  sites. 
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4.0.  Concurrent  Multigrid  Speedup  -  The  Pyramid  Network 


Like  all  local  iterative  algorithms,  the  coupled  depth/slopc  method  takes 
a  long  time  to  converge.  For  realistic  images,  it  literally  takes  thousands  of 
iterations  to  calculate  a  surface.  In  this  section,  we  propose  a  new  concurrent 
multigrid  algorithm  to  speedup  the  coupled  depth/slopc  calculations. 


% 
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Terzopoulos  [1984]  used  multigrid  methods  to  improve  the  convergence  rate 
of  his  surface  reconstruction  algorithm.  This  was  done  for  normal  sequential 
computers.  Multigrid  methods  [Brandt  1977;  Hackbusch  and  TVoltcnbcrg  1981] 
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Figure  -1-8.  Fxninplc  of  fault  tolerance  result 

arc  a  standard  set  of  mathematical  techniques  designed  to  speed  the  convergence 
of  an  iterative  algorithm.  The  basic  idea  is  to  solve  the  problem  Grst  on  a  coarse 
level  and  then  use  that  solution  as  a  first  approximation  for  a  finer  level. 

Multigrid  methods  arc  often  many  orders  of  magnitude  faster  than  their 
single  level  counterparts.  For  instance,  suppose  we  arc  given  an  N  by  /V  grid 
on  which  to  calculate  the  solution  of  a  surface.  The  computation  can  first  be 
mapped  to  a  coarser  grid  of  N/2  by  N/2  mesh  points.  There  arc  two  reasons  for 
the  mulligrid  speedup.  A  single  iteration  on  this  coarse  grid  takes  one-fourth 
the  number  of  the  operations  of  a  single  iteration  on  the  finer  grid.  In  addition, 
since  signals  do  not  have  to  propagate  as  far  on  the  coarse  grid,  it  takes  fewer 
iterations  to  converge  to  a  solution.  So  convergence  on  the  coarse  grid  is  much 
faster  since  fewer  iterations  arc  needed  and  less  computation  is  required  for  each 
iteration. 

When  a  multigrid  algorithm  is  run  on  a  problem  which  has  already  been 
mapped  to  a  parallel  machine,  the  relative  speedup  is  not  as  great  as  a  multigrid 
algorithm  running  on  a  sequential  machine.  This  is  because,  on  a  parallel  ma- 


chine,  the  time  it  takes  to  do  a  single  iteration  on  any  level  is  the  same  no  matter 
which  level  of  coarseness  is  being  solved.  In  the  sequential  case,  iterations  on  the 
next  coarser  level  take  one-quarter  the  speed  of  the  finer  level,  but  in  a  parallel 
machine,  all  the  nodal  computations  are  done  in  parallel  so  all  grids  do  iterations 
equally  fast.  However,  multigrid  methods  still  give  some  algorithmic  speedup  on 
a  parallel  machine  because  it  takes  fewer  iterations  to  converge  on  coarser  levels. 
It’s  just  that  the  overall  relative  speedup  for  introducing  multigrid  methods  is 
not  as  great  as  for  a  sequential  machine. 

There  is  a  problem  with  mapping  these  multigrid  algorithms  to  a  parallel 
computer  like  the  Connection  Machine.  The  multigrid  methods  presently  used 
can  only  process  data  on  one  level  at  a  time.  In  this  sense  these  algorithms 
are  inherently  sequential.  Recently,  Terzopoulos  [1985]  developed  an  algorithm 
which  allows  all  of  the  levels  to  run  in  parallel.  He  does  away  with  the  complex 
coarse-to-fine  and  fine-to-coarse  transfers  of  data  normally  seen  in  multigrid 
algorithms.  However,  Terzopoulos’  algorithm  requires  the  use  of  global  switching 
parameters  which  control  the  information  flow  between  coarse  and  fine  levels. 
During  the  initial  stages  of  the  algorithm,  a  larger  percentage  of  information  will 
flow  from  the  coarser  levels  down  to  the  finer  levels.  As  the  algorithm  progresses, 
more  information  flows  up  from  the  finer  levels,  correcting  the  coarser  levels. 

We  propose  a  better  solution.  We  model  the  mulligrid  problem  as  a  network 
of  resistors  and  capacitors  arranged  in  a  pyramid.  This  pyramid  network  is 
composed  of  meshes  of  resistors  and  capacitors,  where  different  levels  of  meshes 
correspond  to  varying  degrees  of  coarseness.  We  develop  this  new  approach  by 
studying  the  simple  ID  membrane  problem  first.  In  a  membrane  smoothing 
algorithm,  straight  lines  are  fit  to  sparse  depth  constraints  through  a  parallel 
iterative  averaging  process.  Figure  4-9a  shows  sparse  constraints  along  a  line 
with  65  mesh  points.  If  an  initial  guess  of  zeroes  is  assumed,  the  result  after  50 
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a)  sparse  constraints 


b)  after  50  iterations 


c)  after  250  iterations 
Figure  4-9.  Single  grid  example 

iterations,  Figure  4-9b,  is  far  from  convergence.  After  250  iterations,  Figure  4- 9c, 
the  solution  has  still  not  converged.  Convergence  is  reached  after  more  than  500 
iterations.  The  slow  convergence  problem  is  solved  with  the  pyramid  network. 

For  the  ID  membrane  problem,  we  can  model  the  computation  as  a  dis¬ 
tributed  RC  line  as  shown  in  Figure  4-10.  The  exact  values  of  R  and  C  must  be 
carefully  chosen  to  accurately  reflect  the  amount  of  information  that  can  pass 
between  processors  during  each  cycle.  Each  iterative  step  corresponds  to  some 
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small  time  step  for  the  analog  network.  The  effect  of  «a  single  constraint  will  be 
modeled  here  as  a  voltage  step  applied  to  one  end  of  the  line.  The  initial  guess 
is  represented  as  the  initial  voltages  on  all  of  the  capacitors. 


Electrically,  the  line  will  exponentially  converge  to  a  solution  with  some  RC 
time  constant.  Glasscr  and  Dobberpuhl  [1985]  give  an  upper  bound  on  the  time 


constant  as 


r  =  RC 


N{N  +  1) 


with  N  being  the  number  of  nodes  in  the  mesh.  The  speed  of  the  delay  line 
varies  quadratically  with  the  its  length.  An  obvious  speed  up  is  to  add  a  coarser 
RC  line  such  as  that  shown  in  Figure  4-11.  This  line  will  solve  the  problem  four 
times  as  fast  since  it  has  half  the  length  and  the  same  R  and  C  values  as  the 


finer  level.  The  constraint  on  the  finer  level  is  used  as  a  constraint  on  the  coarser 


level  at  the  same  spatial  location. 


Problems  occur  when  constraints  don’t  occur  on  mesh  points,  This  difficulty 
is  resolved  through  careful  selection  of  inter-grid  information  flow.  Suppose  we 
arc  given  constraint  values  at  locations  0  and  7  in  the  nine  node  arrangement 
shown  in  Figure  4-12.  We  use  a  hierarchy  of  four  grids  of  varying  levels  of  reso¬ 
lution  to  solve  the  membrane  problem  here.  We  observe  that  in  this  simple  case, 
wires  connecting  different  levels  need  only  be  unidirectional.  The  direction  of 
information  flow  is  determined  by  the  location  of  the  constraints.  Information 
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Figure  4-11.  Two  level  speedup 


flows  tip  towards  coarser  grids  where  there  are  constraints  and  it  flows  down 
where  there  arc  no  constraints.  This  information  flow  between  levels  is  an  injec¬ 
tion  process.  The  recipient  of  this  information  just  takes  these  values  and  docs 
not  do  any  neighborhood  averaging. 

This  method  of  transferring  constraints  to  coarser  levels  is  very  effective. 
Initially,  only  approximate  constraint  values  are  passed  up  to  the  coarser  grids. 
As  the  problem  is  solved  on  coarser  grids  inform<ation  flows  down  to  the  finer 
grids,  creating  better  constraint  values  to  pass  to  the  coarser  grids.  Tcrzopoulos 
used  a  static  local  weighting  technique  to  pass  constraint  values  to  coarser  grids. 
Ilis  method  has  the  disadvantage  that  the  approximate  constraiut  values  are 
not  improved  over  time  and  so  the  coarser  grids  can  never  have  an  accur.ate 
representation  of  the  surface. 

How  many  coarse  grids  should  be  used?  This  pyramid  algorithm  uses  as 
many  levels  as  it  needs  and  uses  them  ail  in  parallel.  Notice  in  Figure  4-12 
that  the  two  coarsest  levels,  although  built  in,  are  not  used.  Data  only  flows 
into  them  -  no  lower  levels  get  information  from  them.  If  this  surface  patch  was 
adjacent  to  a  sparser  region,  these  coarser  levels  would  then  be  able  to  efficiently 
pass  data  there.  Generally,  the  sparser  the  data  in  a  particular  region,  the  more 
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Figure  4-12.  Intcr-grid  information  flow. 

levels  of  the  pyramid  that  arc  used.  If  data  is  very  dense,  mulligrid  methods 
will  not  enable  any  appreciable  speedup. 

This  pyramid  network  algorithm  has  been  implemented  for  the  ID  mem¬ 
brane  case.  The  solution,  shown  in  Figure  4-13,  took  only  50  iterations.  This 
is  more  than  an  order  of  magnitude  speed  increase  over  the  single  grid  solution. 
The  pyramid  algorithm  converges  even  faster  with  the  addition  of  a  heuristic 
that  uses  the  pyramid  structure  to  develop  a  good  initial  guess. 

We  have  developed  a  concurrent  niultigrid  acceleration  technique  to  speed 
up  the  surface  reconstruction  problem.  This  pyramid  network  is  better  suited 
for  a  parallel  machine  than  conventional  multigrid  techniques  because: 

1.  The  pyramid  network  runs  all  of  the  grids  concurrently. 

2.  The  algorithm  is  much  simpler  to  program  than  standard  multigrid  algo¬ 
rithms. 

The  pyramid  algorithm  has  the  following  advantages  over  Tcrzopoulos’  con¬ 
current  multigrid  approach: 

1.  The  system  does  not  rely  on  any  dynamically  changing  global  parameters. 
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2.  The  coarser  level  solutions  have  the  potential  to  be  more  accurate  since 
the  coarse  constraint  values  improve  gradually  when  the  finer  grid  solution 
improves. 

3.  An  analog  network  implementation  seems  feasible. 

4.  The  system  automatically  chooses  how  many  coarse  levels  .ore  used.  Regions 
where  information  is  sparse  tend  to  use  the  most  grid  levels. 


4.7.  Digital  Alternatives 


Historically,  direct  methods  for  solving  extremely  large  partial  differential 
equation  problems  have  been  prohibitively  expensive.  The  advent  of  locally  con¬ 
nected,  massively  parallel  machines  has  provided  iterative  algorithms  with  even 
more  advantage  over  direct  methods.  Recently,  Pan  and  Reif  [1985]  have  pro¬ 
posed  the  parallel  nested  dissection  algorithm  for  solving  linear  systems.  Their 
algorithm  for  solving  large  sparse  systems  maps  particularly  well  on  a  computer 
such  as  the  Connection  Machine.  Such  a  technique  could  be  applied  to  solve 
the  system  of  coupled  difference  equations  given  by  the  coupled  depth/slope 
method. 

In  general,  direct  methods  are  much  more  difficult  to  implement  and  require 
more  memory  than  iterative  methods.  Iterative  methods  can  easily  make  good 
use  of  approximate  solutions.  This  is  important  in  solving  problems  like  optical 
flow.  The  optical  flow  in  one  time  frame  should  be  very  close  to  the  optical  flow 
in  the  next  frame.  Direct  methods  cannot  make  use  of  good  approximations.  For 
extremely  large  problems,  iterative  methods  are  more  efficient,  that  is,  less  work 
is  required  to  achieve  a  solution  of  a  given  accuracy.  Low  accuracy  solutions, 
can  be  generated  very  rapidly.  The  iterative  multigrid  methods  have  optimal 
time  complexity  rates  for  the  solution  of  linear  systems.  These  multigrid  solvers 
can  also  work  with  the  same  efficiency  on  nonlinear  and  eigenvalue  problems. 
Multigrid  algorithms  have  been  successfully  adapted  to  nonlinear  flow  problems 
such  as  transonic  flows  and  incompressible  Navier-Stokes  equations  [Brandt  81]. 

Terrence  Boult  describes  yet  another  way  to  perform  surface  reconstruction. 
He  uses  continuous  thin  plate  splines  [Boult  1985a!,  but  does  not  attempt  to 
integrate  orientation  information,  does  not  provide  for  discontinuities  of  different 
types  and  makes  no  provisions  for  handling  noisy  measurements.  Therefore, 


from  the  definition  given  in  Chapter  1,  he  does  not  address  the  full  surface 
reconstruction  problem,  only  surface  interpolation. 

Boult  compares  the  gradient  projection  based  algorithm  used  by  Grimson 
with  his  own  method  of  reproducing  kernels  [Boult  1985b].  He  argues  that  the 
reproducing  kernel  method  has  better  parallel  time  complexity  than  the  single 
level  algorithm  of  Grimson.  Assume  k  depth  points  are  given  in  an  m  x  m  mesh. 
Define  the  number  of  grid  points  m2  to  be  n  and  assume  a  typical  10%  sampling 
of  depth  points  so  k  =  n/10.  The  time  consuming  computation  for  this  method 
is  the  solution  of  a  (A:  -  3)  by  (Ik  -  3)  dense  linear  system  of  equations  which 
takes  0(k2)  or  equivalently  0(n2)  time  to  solve  (assuming  zero  communication 
time).  This  compares  favorably  with  O(nlogn),  the  parallel  time  complexity 
that  the  Boult  gives  for  Grimson ’s  algorithm.  The  multigrid  methods  used  by 
Terzopoulos  lowers  this  time  complexity  to  O(n)  on  a  sequential  machine  [Brandt 
1977], 

In  this  analysis,  communication  time  complexity  is  ignored,  despite  the 
fact  that  this  seems  to  be  the  major  bottleneck  in  Boult’s  approach.  The  thin 
plate  spline  generating  algorithm  does  not  exploit  the  locality  that  exists  in 
the  problem.  This  contrasts  with  the  coupled  depth/slope  algorithm.  Surface 
reconstruction  using  iterative  multigrid  methods  has  better  computation  and 
communication  time  complexities. 

Recently  Jim  Clark  has  developed  a  novel  approach  to  surface  interpola¬ 
tion  [Clark  1985].  His  idea  consists  of  three  steps.  First,  he  maps  the  sparse 
sensory  data  to  regularly  spaced  points  in  the  2D  grid.  Second,  he  uses  stan¬ 
dard  techniques  from  sampling  theorems  in  digital  signal  processing  and  does 
a  single  convolution  of  regularly  spaced  points  to  recreate  a  surface.  Finally, 
an  inverse  mapping  is  performed  to  map  the  reconstructed  surface  back  to  the 
original  domain.  Unfortunately,  his  work  does  not  extend  to  handle  arbitrary 


discontinuities.  Furthermore,  his  method  fails  to  extrapolate  the  unconstrained 
edges  of  the  surface  as,  for  example,  the  biharmonic  operator  would. 

4.8.  Summary 

This  chapter  has  developed  a  digital  algorithm  based  upon  the  coupled 
slope  depth  model  to  solve  the  surface  reconstruction  problem.  The  algorithm  is 
extremely  well  suited  for  a  locally-connected  massively-parallel  computer  such 
as  the  Connection  Machine.  Since  this  algorithm  comes  directly  from  the  cou¬ 
pled  depth/slope  model,  the  algorithm  can  handle  constraints  and  discontinu¬ 
ities  at  any  order  of  derivative  of  the  surface.  In  addition,  the  digital  coupled 
depth/slope  algorithm  was  shown  to  have  the  following  advantages  for  a  Con¬ 
nection  Machine  type  implementation: 

1.  It  requires  only  nearest  neighbor  communication  between  processors.  No 
information  from  processors  other  than  the  four  nearest  neighbors  is  ever 
needed.  This  includes  depth  and  orientation  constraints,  discontinuity  lo¬ 
cations  as  well  as  current  parameter  values. 

2.  The  algorithm  was  shown  to  be  extremely  fault  tolerant.  Broken  processors 
merely  cause  a  loss  of  resolution;  the  solution  is  still  roughly  correct. 

3.  Each  node  requires  only  fixed  point  additions,  subtractions  and  averages. 
This  constraint  is  extremely  important  for  a  computer  such  as  the  Connec¬ 
tion  Machine  which  uses  simple  one-bit  processors. 
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CHAPTER  5 


ANALOG  IMPLEMENTATION 


The  coupled  depth/slope  algorithm  discussed  in  Chapter  4  runs  much  faster 
when  implemented  with  analog  circuitry.  Instead  of  simulating  an  analog  net¬ 
work  with  a  digital  computer,  this  chapter  describes  possible  analog  implemen¬ 
tations  of  the  coupled  depth/slope  network.  First,  we  discuss  how  to  build 
the  necessary  subtractor  constraint  device  and  then  we  discuss  how  to  set  the 
constraint  values. 

5.1.  Building  the  Subtractor  Constraint  Element 

The  hardest  part  of  building  the  coupled  depth/slope  network  is  the  con¬ 
struction  of  the  subtractor  elements  which  are  responsible  for  the  coupling  be¬ 
tween  the  slope  and  depth  planes.  As  Chapter  3  states,  any  of  the  terminals  of 
this  device  can  be  inputs  or  outputs.  One  way  of  building  this  device  is  with 
an  ideal  transformer  in  the  configuration  shown  in  Figure  5-1.  Here  we  rely  on 
the  mutual  inductance  of  the  two  coils  to  produce  the  desired  coupling  effect 
[Senturia  and  Wedlock  1975,  p.  135].  With  a  turns  ratio  of  one,  the  voltages 
across  the  two  inductors  will  be  the  same,  so  A  -  B  =  C.  The  transformer 
circuit  requires  alternating  current  to  function  and  dissipates  no  power. 

The  use  of  transformers  makes  the  coupled  network  equivalent  to  the  elastic 
plate  analogy  built  by  MacNeal  [1951]  in  the  early  1950’s.  His  elastic  plate 
analogy  network  contained  twelve  nodes  and  could  be  extended  to  solve  problems 
such  as  beam  deflections  and  transient  vibrations.  There  are  several  problems 
with  any  transformer-based  network.  The  most  severe  problem  is  building  error- 
free  transformers.  A  second  problem  is  that  it  is  virtually  impossible  to  build 


a  reasonable  transformer  on  a  VLSI  chip.  This  precludes  the  possibility  of 
constructing  a  large  scale  system. 


Another  way  of  constructing  the  ideal  sublractor  element  is  with  the  help 
of  negative  resistors  as  shown  in  Figure  5-2.  By  Kirchoff’s  current  law  the  sum 
of  the  currents  into  the  central  node  node  X  is  sero,  so 
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or  rearranging  terms. 


A-  B  =  C 


There  are  a  number  of  ways  of  building  negative  resistances.  One  possibility 
b  to  use  a  careful  combination  of  capacitors  and  inductors  as  shown  in  Figure 
5-3.  FYom  above,  we  require: 


Zl  =  -Zc 
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Figure  5-2.  Negative  resistance  implementation. 


Figure  5-3.  FYcqtiency  domain  implementation. 

Not  only  is  alternating  current  required,  but  the  frequency  has  to  exactly  match 
the  values  of  L  and  C.  This  is  a  difficult  circuit  to  build  correctly,  and  again  we 
cannot  build  it  onto  an  integrated  circuit. 

Another  implementation  strategy  is  to  implement  the  device  with  two  ana- 


log  subtractors  and  one  analog  adder,  as  was  shown  in  Figure  3-2.  An  adder  is 
needed  to  compute  B  +  C  and  the  two  subtractors  are  used  to  compute  A  -  C 
and  A  -  B.  Adders  and  subtractors  are  simple  to  build  with  operational  am¬ 
plifiers.  The  full  subtract  constraint  element  circuit  is  shown  in  Figure  5-4.  A 
working  version  of  this  circuit  has  been  built  and  tested. 

Problems  occur  in  real  life  when  we  try  to  construct  ideal  devices  out  of  non¬ 
ideal  parts.  Each  op  amp  adder  and  subtractor  circuit  has  tiny  errors  in  resis¬ 
tance  values,  finite  gain,  finite  input  impedance  and  non-zero  output  impedance. 
These  errors  cause  slight  inaccuracies  in  the  calculated  values  and  the  constraint 
A  —  B  =  C  will  not  be  true  in  general.  The  small  error  doesn’t  bother  us  from 
an  algorithmic  standpoint  since  the  inaccuracies  are  very  small  and  are  averaged 
together  in  each  constraint  element.  There  is  a  problem  from  the  electrical  point 
of  view,  though.  The  circuit  described  sets  the  output  of  all  of  its  op  amps  to 
values  which  are  slightly  different  from  their  normal  output.  Since  op  amps  have 
an  extremely  small  output  impedance  a  lot  of  current  will  flow,  possibly  burning 
the  op  amp.  The  problem  is  reduced  by  explicitly  increasing  the  op  amp  output 
impedance  with  an  extra  resistor  in  the  feedback  path. 

Perhaps  a  better  approach  is  to  implement  directly  the  negative  resistance 
values  given  in  Figure  5-3  with  negative  impedance  converters  [Siebert  1986], 
With  two  negative  impedance  circuits,  the  subtract  constraint  element  can  be 
constructed  with  two  op  amps  and  eight  resistors  as  shown  in  Figure  5-5.  Also, 
the  rivalry  and  error  problems  described  for  the  circuit  in  Figure  5-4  aren’t 
present  in  this  implementation. 

The  big  advantage  of  these  op  amp  circuits  is  that  operational  amplifiers 
are  readily  fabricated  with  VLSI  technology.  The  idea  of  surface  reconstruction 
on  a  chip  is  not  too  far  in  the  future.  The  op  amps  need  not  consume  excessive 
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Figure  5-4.  Operational  amplifier  implementation. 


chip  space  since  the  the  DC  biasing  circuitry  can  be  shared  among  many  op 
amps. 

5.2.  Setting  Constraint  Values 

In  Chapter  3,  we  proposed  to  set  constraint  values  with  a  voltage  source 
connected  through  a  resistor  to  the  constraint  mesh  point,  the  Thevenin  equiv¬ 
alent,  see  Figure  5-6a.  The  depth  measurement  u  corresponds  to  the  value  of 
the  voltage  source  and  the  conductance  of  the  resistor  is  a,  with  a  defined  as 
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Figure  5-5.  Negative  impedance  converter  implementation. 


a  measure  of  the  confidence  in  the  measurements.  The  more  typical  connection 
proposed  in  surface  approximation  networks  [Marroquin  1985;  Koch,  Marroquin 
and  Yuille  1986;  Poggio  and  Koch  1985;  Hutchinson  1986]  is  to  use  the  Norton 
equivalent  circuit,  shown  in  Figure  5-6b. 


The  Thcvcnin  equivalent  used  to  set  constraints  has  several  important  im¬ 
plementation  advantages. 


1.  While  resistor  meshes  .are  generally  insensitive  to  resistor  value  variations 
[Karplus  1958],  the  Norton  circuit  is  very  sensitive  to  the  resistance  varia¬ 
tions  that  are  common  over  VLSI  process  corners  [Hutchinson  1986].  The 
voltage  seen  at  u  varies  directly  with  variations  in  resistor  values.  The  sim¬ 
ilar  problem  docs  not  occur  for  the  Thevenin  circuit.  Variation  in  a  resistor 
value  simply  changes  the  confidence  in  u. 


2.  Using  the  Thevenin  equivalent  circuit  also  has  the  advantage  of  reduced 

power  consumption.  We  want  to  maximize  the  voltage  swing  in  the  resistor 
mesh  to  achieve  the  largest  resolution  possible.  Power  dissipation  is  a  prob¬ 
lem,  therefore,  since  power  is  We  set  the  resistor  values  in  the 

mesh  as  large  as  possible  to  minimize  circuit  power.  It  can  be  shown  that 
the  Norton  circuit  dissipates  much  more  power  than  the  Thevenin  circuit 
when  driving  a  large  resistive  load. 

3.  Finally,  it  is  easier  to  implement  the  Thevenin  circuit  since  the  variable  a 
parameter  sets  only  the  resistance  value  and  is  not  involved  with  the  value 
of  the  voltage  source.  The  a  value  plays  a  part  in  both  the  current  source 
and  the  resistance  values  for  the  Norton  equivalent  circuit. 

Unfortunately,  it  is  simpler  to  build  current  sources  than  voltage  sources 
with  today’s  silicon  technology.  Nonetheless,  the  Thevenin  circuit  has  the  po¬ 
tential  for  more  process  variation  immunity,  less  power  dissipation,  and  simpler 
overall  system  implementation. 

5.3.  Alternative  Analog  Networks 

There  are  other  possible  analog  networks  which  implement  biharmonic  smooth- 


ness.  The  biharmonic  mask  can  be  synthesized  electrically  by  using  the  tricks 
discussed  in  the  last  section  for  manufacturing  negative  resistance  components. 
For  example,  the  biharmonic  mask  can  be  realized  in  the  frequency  domain  with 
a  network  of  resistors,  capacitors  and  inductors  [Volynskii  65].  Figure  5-7  shows 
a  typical  circuit.  Unfortunately,  these  networks  are  hard  to  build  in  today’s 
VLSI  technology  because  inductors  are  difficult  to  fabricate.  Also,  small  errors 
in  capacitance  and  inductance  values  will  cause  large  variations  in  the  solution. 
A  more  plausible  network  would  take  advantage  of  the  negative  impedance  con¬ 
verters  used  in  Figure  5-5. 

The  dual-Poisson  network  described  in  Chapter  3  can  be  implemented 
straightforwardly.  However,  depth  values  (or  slope  values)  must  be  given  at 
each  point  around  the  entire  boundary.  This  network  cannot  implement  the 
biharmonic  frcc-edge  condition  required  for  surface  reconstruction.  Also,  there 
can  be  no  discontinuities  in  the  interior  of  the  image  unless  the  slope  or  depth 
is  known  at  every  point  on  both  sides  of  all  discontinuities.  The  dual  Poisson 
network  may  be  .appropriate  for  restricted  domains  but  is  not  useful  for  general 


surface  reconstruction  applications. 


5.4.  Summary 

We  have  proposed  implementation  strategies  for  building  an  analog  ver¬ 
sion  of  the  coupled  depth/slope  method.  The  rough  estimated  speed  of  such  a 
network  is  on  the  order  of  microseconds.  Unfortunately,  this  network  will  be 
plagued  by  the  problems  common  to  all  analog  networks:  I/O  boundedness, 
inflexibility  and  limited  resolution.  Clearly,  the  speed  of  the  system  will  be  de¬ 
termined  by  how  fast  data  can  be  loaded  and  unloaded  into  the  network.  It  is 
very  difficult  to  modify  the  algorithms  when  they  are  comped  into  hardware. 
Noise  limits  the  precision  of  voltages  stored  on  nodes  in  the  system.  Because  of 
these  problems,  analog  networks  have  been  little  more  than  a  curiosity  in  the 


CHAPTER  6 

FUTURE  DIRECTIONS 
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The  coupled  depth/slope  model  provides  a  new  way  to  think  about  surface 
reconstruction.  As  always,  new  solutions  tend  to  open  new  areas  of  research. 
Previously,  we  ignored  the  black  boxes  that  performed  discontinuity  detection, 
“shape  from”,  and  stereo.  We  assumed  that  the  results  of  these  early  vision 
procedures  were  available  as  input  to  the  surface  reconstruction  module.  An  area 
for  future  research  is  to  combine  these  early  vision  processes  with  the  surface 
reconstruction  module  to  improve  the  performance  of  the  overall  system. 

6.1.  Combining  Vision  Modules 
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Discontinuity  detection  is  an  active  area  of  research  in  computer  vision  to¬ 
day.  Terzopoulos  [1985]  and  Marroquin  [  1985]  both  looked  at  combining  discon¬ 
tinuity  detection  with  surface  reconstruction.  Terzopoulos  developed  a  heuristic 
for  recognizing  discontinuities.  After  a  surface  has  been  fully  reconstructed,  he 
looks  for  significant  changes  in  curvature  in  order  to  establish  the  placement 
of  discontinuities.  Marroquin  used  binary  line  processes  to  represent  surface 
discontinuities.  He  minimizes  a  non-quadratic  cost  function  to  locate  disconti¬ 
nuities. 

Either  of  the  above  methods  can  be  incorporated  into  the  coupled  depth/slope 
model.  Terzopoulos’  method  can  be  implemented  by  building  the  resistors  in 
the  resistor  mesh  as  fuses.  The  resistors  open  circuit  when  the  current  passing 
through  them  exceeds  some  threshold.  This  threshold  should  be  set  to  some  fi¬ 
nite  value  only  at  the  locations  of  edges  in  the  primal  sketch.  This  will  decrease 


the  possibility  of  false  discontinuity  detections.  In  this  way,  discontinuities  at  any 
level  of  smoothness  could  be  detected.  With  the  concurrent  multigrid  approach, 
discontinuities  at  different  levels  of  resolution  can  be  detected  and  represented. 

Similarly,  the  stereo  module  could  be  effectively  combined  with  the  coupled 
depth/slope  algorithm.  Stereo  modules  typically  use  a  hierarchical  style  ap¬ 
proach  to  solve  the  stereo  matching  problem.  The  simpler,  more  obvious  matches 
are  done  on  coarser  grids.  These  coarser  matches  are  then  carried  through  onto 
the  finer  grids.  If  the  surface  reconstruction  module  were  integrated  and  per¬ 
formed  concurrently  with  the  stereo  module,  the  following  advantages  would 
be  seen.  First,  coarser  surface  reconstruction  could  start  as  soon  as  the  stereo 
module  registered  some  coarse  matches  with  disparities.  This  would  speed  the 
running  time  of  the  surface  reconstruction  module.  Second,  coarse  reconstructed 
surfaces  could  be  fed  back  to  the  stereo  module  to  aid  matching  at  finer  levels. 
This  will  help  the  stereo  module  to  find  more  depth  points  with  fewer  errors, 
resulting  in  a  faster  convergence  and  a  more  accurate  reconstructed  surface. 

6.2.  Other  Applications 

The  coupled  depth/slope  network  is  not  restricted  to  solving  only  the  surface 
reconstruction  problem.  This  network  generalizes  to  arbitrary  dimensions  and 
arbitrary  levels  of  smoothness.  For  example,  we  have  seen  that  Laplacian,  Pois¬ 
son,  and  biharmonic  type  problems  are  easily  solved  with  this  network.  These 
partial  differential  equations  arise  in  such  computationally  demanding  fields  as 
hydrodynamics,  aerodynamics,  weather  forecasting  and  structural  analysis.  The 
hardware  described  in  Chapter  5  is  being  studied  for  use  in  thin  beams  or  can¬ 
tilever  problems  in  structural  analysis.  A  simple  example  is  shown  in  Figure  6-1 
where  an  elastic  beam  is  fixed  on  its  right  and  left  edges  and  a  downward  force 


Figure  6-1.  Elnstic  Beam  Problem 


is  applied  at  the  center  of  the  beam.  The  final  displacement  of  the  beam  must 
be  found. 

The  coupled  dcplh/slopc  network  Ciisily  handles  problems  such  as  these. 
The  input  forces  arc  fed  as  input  at  the  bottom  of  the  network  (sec  Figure 
6-2)  using  appropriately  valued  voltage  sources.  There  is  no  need  for  any  of 
the  resistor  meshes  because  these  are  dense  constraints.  No  minimization  is 
necessary;  only  integration  is  performed.  Extra  constraints  are  necessary  to 
regain  the  constants  of  integration.  Alternatively,  we  could  feed  dense  deflection 
data  into  the  top  of  the  array  and  read  the  forces  from  the  bottom.  Or  we  could 
feed  only  sparse  data  .and  add  a  resistor  mesh  for  some  appropriate  level  of 
smoothness.  In  summary,  we  can  give  the  network  any  data  that  we  have  about 
the  moments,  forces,  slopes,  or  deflections  of  the  beam.  These  constraint  values, 
with  their  corresponding  confidence  values,  .are  incorporated  into  the  final  beam 
solution.  We  can  directly  read  dense  solutions  of  the  loading,  the  shearing  force, 
the  bending  moment  and  the  beam  deflection. 


0.3.  Problems  to  be  addressed 


The  coupled  dcpth/slope  network  has  raised  a  number  of  questions  which 
remain  to  be  answered. 


Figure  6-2.  Elastic  Beam  Solution:  a)  Loading  diagram  b)  Shearing  force  diagram 
c)  Bending  moment  d)  Deflection  liue 

IIow  similar  is  the  post-smoothing  approach  (Section  4.4)  to  the  conven¬ 
tional  least  squares  averaging  approach? 

This  report  has  concentrated  on  parallel  implementations  of  the  coupled 
depth/slopc  network.  How  well  does  this  algorithm  work  with  a  Gauss 
Seidel  iterative  scheme  on  a  sequential  scheme?  Does  the  method  have 
a  low  enough  smoothing  rate  to  be  appropriate  for  a  standard  multigrid 
implementation? 

The  concurrent  multigrid  algorithm  described  in  Section  4.0  has  been  im¬ 
plemented  for  solving  one-dimensional  membrane  problems.  It  must  be 
modified  for  the  2D  membrane  case  and  for  the  coupled  depth/slopc  algo¬ 
rithm.  What  are  the  exact  speeds  of  these  algorithm  on  the  Connection 
Machine  and  on  a  pyramid  machine?  How  do  these  rates  compare  to  the 
speedup  afforded  by  conventional  multigrid  algorithms? 


The  2D  form  of  the  generalized  coupled  depth/slope  network  must  be  de- 


veloped.  It  should  be  possible  to  show  convergence  for  both  the  analog 
network  and  its  digital  implementation. 


What  is  the  best  way  to  build  the  analog  subtract  constraint  device?  Will 
the  circuits  proposed  in  Chapter  5  converge  to  a  solution  or  oscillate  and 
diverge? 


CHAPTER  7 


CONCLUSION 


In  summary,  we  have  developed  the  coupled  depth/slope  model  for  surface 
reconstruction.  This  is  a  computational  model  with  a  number  of  advantages  over 
earlier  models.  If  only  depths  are  given,  the  model  reconstructs  a  surface  with  the 
same  biharmonic  smoothness  used  by  Terzopoulos  and  Grimson.  In  addition,  as 
depths  are  computed,  orientation  constraints  can  be  incorporated  with  no  extra 
work.  Discontinuities  in  depth  and  orientation  can  be  used  to  compute  the 
reconstructed  surface.  Furthermore,  the  coupled  depth/slope  model  generalizes 
to  handle  constraints  and  discontinuities  of  any  order  derivative. 

The  model  appears  to  be  biologically  feasible  in  that  it  fits  the  four  cri¬ 
terion  given  in  Chapter  1:  parallelism,  uniformity,  locality  and  fault  tolerance. 
Its  biharmonic  smoothness  properties  are  observed  in  psychophysical  tests  of 
humans.  Also,  Marr  argues  that  the  human  2^-D  sketch  redundantly  represents 
both  slopes  and  depths.  This  compares  favorably  with  the  dense  depth  and  slope 
representations  generated  by  the  coupled  depth/slope  network. 

Besides  appearing  to  be  biologically  feasible,  the  coupled  depth /slope  model 
leads  naturally  to  compuational  solutions  on  both  digital  computers  and  in  ana¬ 
log  circuitry.  In  fact,  a  digital  algorithm  based  upon  the  coupled  depth/slope 
model  was  implemented  and  was  found  to  be  extremely  well  suited  for  a  locally- 
connected  massively-paralle!  computer  such  as  the  Connection  Machine.  Since 
this  algorithm  comes  directly  from  the  coupled  model,  the  algorithm  can  handle 
constraints  and  discontinuities  at  any  order  of  derivative. 

The  coupled  depth/slope  algorithm  was  shown  to  have  the  following  ad¬ 
vantages  for  a  Connection  Machine  implementation.  It  required  only  nearest 
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neighbor  communication  between  processors.  No  information  from  processors 
other  than  the  four  nearest  neighbors  is  ever  needed.  This  includes  depth  and 
orientation  constraints,  and  discontinuity  locations  as  well  as  current  depth  and 
slope  values.  The  algorithm  was  shown  to  exhibit  a  high  degree  of  fault  toler¬ 
ance.  Broken  processors  cause  only  a  loss  of  resolution;  the  overall  solution  is 
still  correct.  Each  node  requires  only  fixed  point  additions,  subtractions  and 
averages.  This  constraint  is  extremely  important  for  a  computer  such  as  the 
Connection  Machine  which  uses  simple  1-bit  processors. 

In  addition  to  the  above  results,  we  have  suggested  a  new  concurrent  multi¬ 
grid  strategy  for  speeding  up  a  massively  parallel  digital  implementation.  The 
pyramid  network  allows  iterations  to  be  performed  on  every  level  simultaneously. 
The  implementation  is  much  simpler  than  that  of  the  conventional  multigrid 
algorithm.  The  system  converges  to  a  solution  without  the  addition  of  any  dy¬ 
namically  changing  parameters.  The  coarser  level  solutions  are  potentially  more 
accurate  since  the  coarser  constraint  values  improve  gradually  when  the  finer 
grid  solution  improves.  An  analog  implementation  of  the  pyramid  network  is 
plausible. 

In  terms  of  biological  vision  systems,  our  computational  model  for  surface 
reconstruction  possesses  remarkable  similarities  to  the  human  vision  system. 
It  remains  to  be  seen  if  biological  systems  actually  do  create  a  dense  surface 
description  or  whether  they  perform  some  more  complicated  process. 

In  terms  of  the  construction  of  artificial  machine  vision  systems,  it  is  quite 
clear  that  we  will  be  able  to  build  systems  which  can  perform  surface  recon¬ 
struction  at  video  scan  rates.  Unfortunately,  no  artificial  later  vision  algorithms 
exist  that  can  make  good  use  of  this  immense  volume  of  data  at  these  fast  rates. 
In  effect,  we  have  shifted  the  vision  bottleneck  from  early  to  immediate  vision 
and  we  are  now  pushing  it  further  towards  later  vision. 
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